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AN ANALYTICAL AND EXPERIMENTAL STUDY OF THE INFLUENCE 
OF COMPONENT VISCOELASTIC SHEAR MODES ON 
LAUNCH VEHICLE TORSIONAL DYNAMICS 
’By 

Ross L. Goble 
ABSTRACT ' 

* This research report presents' an analytical investigation ©f the 
influence of viscoelastic shear modes of a solid propellant motor- on the 
torsional -dynamics of the parent launch vehicle ‘system. The equations 

, * v * 

of motion which describe the behavior of the propellant in pure shear 
are derived from the theory of elasticity/ with the propellant 
mathematically modeled as a thick-walled' hollow cylinder. The resulting 
general solution of the governing Bessel’s equation is then evaluated 
for free-free,, free-pinned, and pinned-free boundary conditions . The 
frequency determinant is also established for the same boundary condi- 
tions. A correspondence principle is Invoked to obtain the associated 
propellant frequency characteristics in the viscoelastic regime. -The . 
composite system problem consisting of the combined propellant and 
launch vehicle components is then developed using Lagrange's equations, 
for dissipative systems with coupling constraint conditions expressed as 
Lagrange Multiplier relations. Small damping assumptions are made, and 
the problem- is expressed in matrix form.. A frequency scaling technique 
is used to circumvent the singularity problems which arise in the 
solution of the system eigenvalue problem due to the presence of the 



Lagrange multiplier relations. A simplified experimental technique is 
then employed to verify the method of systems analysis. 

Application of the analysis to the torsional dynamics of a con- 
temporary axisymmetric launch vehicle with a small solid propellant 
segment attached as an upper stage is presented. Comparisons of system 
dynamic characteristics for the assumptions of elastic, viscoelastic 
and rigid propellant inertias are provided. Results indicate the 
necessity for considering propellant shear flexibility about the 
rotational axis of symmetry in assessing composite system dynamic 
response to oscillatory torsional excitations during flight. 
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I. INTRODUCTION 


Various analytical methods have "been developed for examining the 

dynamical behavior of viscoelastic propellants to include the effects 

* 

of temperature and pressure [l-^] , angular oscillation effects on 
ablation [5 ] , thermomechanical behavior, during cyclic loading [6j , free 
and forced njultiaxis vibrations jj-lS] , and other dynamic phenomena 
[l 9 .j 26] . In short, the viscoelastic cylinder, case' bonded at its outer 

* ’ > i 

boundary to an elastic cylinder.) has been rigorously treated analyti- 

f ’ 1 

cally for a. wide range of significant boundary conditions and. motion 

1 * s- » 

characteristics. These analytical methods deal in general with highly 
idealized structures, but have application to ‘launch vehicles wherein 
the solid propellant motor is the prime mover. Even so, Baltrukronis 
[21]. has pointed out the lack of well-founded methods to quantitatively 
evaluate the contributions of the propellant to the dynamic responses 
of the composite structure,. Achenbach [22] also has stated that the 
interaction of the combustion process and the mechanical vibrations of 
the solid rocket should be further investigated. 

An equally important interaction is that arising in launch vehicle 
systems where small solid propellant rocket motors are used in upper 
stages of launch vehicles in which the main stage is liquid-fueled. In' 
this case the solid propellant prior to its ignition is generally 
treated as a rigid body in analyzing the overall vehicle modal 
characteristics. Otherwise, the coupled system thus described has 
received little attention in the literature. The rigid body simulation 
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of the solid propellant , however, is not thought to he sufficient for - 
analysis of the vehicle system response to< oscillatory torsional 
excitations which can occur in the flight profile. The purpose of the 
present paper is to investigate the torsional response characteristics 
of a more realistic launch vehicle system comprised of a long elastic 
thin-walled vehicle cylinder containing a concentric viscoelastic solid 

x s 

propellant upper stage attached over a short portion of the - elastic 
vehicle cylinder. Comparisons of the response frequencies , and mode 
shapes of the total system for various assumptions concerning the visco- 

* T 

elastic cylinder are presented. The primary value of the investigation 
lies in the fact that the majority of the work done- 1 to date involving 
elastic and viscoelastic cylinders with application to solid propellants 
has pertai-ned to the solid propellant and its bonded container per se 
without relation to the overall systems problem. 

As pointed out by Achenbach [23]' , analytical methods for determina- 
tion of prope.llant dynamics are important, not. only to ensure propellant 
structural integrity, but also to provide information on propellant 
natural frequencies and mode shapes for analysis of the complete launch 
vehicle dynamics. It is a problem in this latter category to which the 
present investigation is addressed. 

As previously mentioned, while the dynamics of solid propellants 
have been examined considerably in the literature, natural vibration 
characteristics of propellants in twist or torsion appear to have 
received the least attention. 
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There have been some recent works, however, which treat the torsion 

I 

problem well* Ail example Is the work by Kelkar which provides a 

* 4 

closed-form solution for a hollow elastic, cylinder bonded to a thin 

* - * 

casing of a different material. Three-dimensional elasticity theory is 

» 

■utilized, and the approach represents a higher degree of sophistication 

* 

than torsional solutions such as Achenbach-' s • [25] which involved only 
one displacement component. 

Other investigators such as Sarm and Shaffer' £26] and Baltrukronis, 
Chi, and Gottehberg £27] performed elastic transverse vibration studies 
for layered cylinders in which two displacement components were obtained. 
Since the problem under consideration, however, is one of total vehicle 
torsional response in which the encased' propellant represents but one 
entity, it was important to survey the literature for studies pertinent 
to dynamics of vehicle systems,. It was found that even when launch 
vehicle systems have been considered, attention has generally been 
directed analytically and experimentally toward modes of vibration and 
response characteristics which do not pertain to torsion [28-53] « 

Matrix methods for discrete mass launch vehicle representation are, of 
course, plentiful and practical since such vehicles are .generally 
composite axisymmetric bodies, discontinuous with regard to distribution 
of physical characteristics along the axis of symmetry. (An example is 
provided in reference [3^ • ) 'While many analytical methods exist for 
motion analysis of such vehicles, using single beam discrete mass 
representation [35] as well as branched-beam methods [36J , recourse- to 
discrete analogies and numerical methods of solution via computer are 
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commonplace; therefore, existing analogies will be used for the dis- 
continuous elastic cylinder portion of the system investigated herein. 

Governing equations' :for' the pure shear natural modes and 
frequencies of a thick-walled’ cylindrical prqpellant are derived from 
basic linear elasticity theory considering free-free, free-pinned, and 
pinned-free boundary conditions. A dynamical correspondence’ principle 

J J ’ 

1 J ' 

is then used to obtain the natural frequencies in the viscoelastic 
medium. The short viscoelastic cylindrical propellant - long elastic 
cylindrical vehicle system' is coupled through the use of Lagrange 
Multiplier constraint relations. Lagrange *'s equations .are then used to 
develop the total system eigenvalue problem. Variations in the elastic 
restraints at the viscoelastic disc - elastic cylinder interface then 
permit system dynamic characteristics to be examined on a parametric 
constraint condition basis., 

The eigenvectors obtained define the modal coefficients of the 
component modes which are subsequently used to establish system modes • 
at the composite system frequencies. The analytical method thus 
developed is sufficiently general to handle a wide array of linear non- 
homegeneous systems of equations wherein the physical dynamics problem 
includes elastic or viscoelastic elements separately or in combination. 
As previously mentioned, the method has particular application to the 
launch vehicle system problem at hand, 

A simplified, disc -shaft -experimental model used to verify the 
•method of systems analysis is discussed, and analytical system natural 
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frequency characteristics are compared against the experimentally 
determined data. 

The method of analysis is then applied to a characteristic launch 
vehicle. Elastic, viscoelastic , and rigid propellant shear modes are used, 
respectively^ in combination with launch vehicle torsional modes to pro- 
vide a comparative study as to the effects of such propellant 

. v « 

assumptions on launch vehicle system torsional response characteristics. 



II. DERIVATION OF THEORY 


The approach taken in solving the systems dynamical problem is to 
first develop the dynamic elasticity solution for the propellant segment 
(figure 1) in a pure shear mode. (Justification for the assumption that 
the propellant segment behaves in pure shear is provided later in the 
-discussion of the theory application to a launch vehicle, system 
problem.) * 

Both the elastic and viscoelastic propellant dynamic characteristics 
are developed in the present chapter. 


Pure Shear Elasticity Solution 

The tangential equilibrium equation for plane problems in polar ' 

’ , 1 * 1 ^ f 

coordinates is [37]* . . , ' - 


1 + d T rcp 

r dcp dr 

The body force per unit volume, I)p, 
problem as 


+ 2 = 0 . - (2-1) 

a > 

» *■ 

can be- expressed' for the shear 




d^p 

dt^ 


( 2 - 2 ) 


but p = cpr, thus 

Sq> = (2-3) 

Now, making the assumption of pure shear motion, substituting relation 
(2-3-)- into (2-1), and simplifying, the governing equation becomes 


6 






8 


OT r 


5r 2 qp - 2 t^ - r -gj* = 0 


Expressing the shear stress in terms of the shear modulus , 


(2-JO 


^Vcp ~ Of? rep 


and it can be seen from figure 1 that the shearing strain is 


, > 3cp 

?rqr = r 5? 


(2-5) 


( 2 - 6 ) 


Therefore 


Y rcp = % r S? 


(2-7) 


Also 


St- 


rep 


Sr 




( 2 - 8 ) 


Making the substitution of relations ( 2.-7 ) arid .(.2-8), into (2-'4) and 
simplifying, the governing equation takes the form 

= ° («) 

Now p .= cpr, and if it is assumed that the cylinder oscillates with 
undamped harmonic motion, p(r,t) = p(r)e^ 0 ^, and the first order 
derivative of equation (2-9) is obtained as 


a 


8cp o r 3r ~ p _ 1 8p _ _P_ = e icnt/Pl _ 


Sr 


5? * ^ 


V-i) 


(2-10) 


where for the sake of brevity, p = p(r,t) and p = p(r) here and in 
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the following development. The second order derivative of (2-9) is 
developed from (2-10) as 


d 2 p dp 


„2 dp ' 


0. . i_^ + _.Sl_ ^ + |) (M1) 


Also 


T 


= £ B -O^e^ £ 


( 2 - 12 ) 


Making substitution of relations (2-10 ) , (2-il ) } and (2-12) into 
equation (2-9 -) f and multiplying through by results in the 

differential equation 


P" P* P 1 | 2P | 3P* 3P , _5_ Q 2 P _ 0 

r^ r^ r^ r 


r r 2 ^ 


(2-13) 


which can be simplified to 


p" + El + 
* r 


(D 


2 8 


% r 2 


P = 0 


Mating the substitution 'E-^- = A 2 into (2-1^) yields 

g f 


+ e: + (# . i) P . o 


(2-14) 


(2-15) 


which is precisely the Bessel’s equation of first order with parameter 
A, for which the solution is 


p(r) = AJ 1 (Ar) + BY 1 (Ar) 


( 2 - 16 ) 
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The physical boundary conditions for the free-free natural mode 
investigation are that the shear vanish at r - R 0 and r = Rjj 


therefore. 


Sep i- _ Sep 


Srj R , Sr 


- 0 (by equation ('2-7),) j hence — = - 

Eo dr r 


(by equation (2-10) )j thus 


dp 

dr' 


p(Ro) dp- p(Ri) 

* 3r 

“i 


Ro *o ' U * Ri R i 


(2-17) 


Substituting equation (2-l6) into the outer and inner radius boundary •- 
conditions of (2-17) yields - * 

p(Ro) 


AJ-l’CARq) + BYx^ARo) = 


ARo 


AJl'(ARi) + BY]_‘ (?®i) = 


P(%) 

ARj_ 


(2-18) 


(2-19) 


Substituting for p(Ro) and p(R-j_) in equations (2-18) and (2-19) 

terms of the exact solution (2-16) and rearranging, the resulting set of 

, ^ / 
simultaneous equations is : 

1 


aJj X '(ARo) - J^ARojj + B^'CARo) Yq^) 

A^^ARi) - ~ J!(ARi)j + B^iHARi) - Y 1 (AR i 


= 0 - • 


= 0 ( 2 - 20 ) 


Now the recurrence relations for Bessel functions include the following: 

r y J 

* ’ 

xJ n '(x) =xJ n . x (x) - nj n .(x) 


where -n indicates the order, and x, the argument. Also 


( 2 - 21 ) 
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xl n '(x> = xY n _ L (x)‘ - nY n (x) 


(. 2 - 22 ) 


Making substitution of (2-21) and ('2-22) into (2-20)', the resulting 

. ' - 1 

set of equations, simplified and expressed in matrix form is /* 


AR o J 0 (AE o ) - 2J i; (ARo) 7<R o y 0 CKR o ) - ^(ARq) 

e 

AR i J 0 (AR i ) - 2J 1 (AR i ) ^iY 0 (ARi) - 2Y 1 (AR i ) 


r, \ 

A 


>- 


o . 


B 

V J 


(2-25) 


The system of equations (2-25) has a solution other than the 
trivial solution if, and only if, the determinant of the coefficients 
is equal to zero. The resulting secular equation for the roots A n is 
thus 


f 


An^JotARb) - 2J 1 (AE Q ) 
RoYotARo) - 2Y 1 (AE 0 ) 


A n : 


)j ^A n R i Y 0 (AR i ). - 2Y 1 (AR i )] 

)j ^A n RiJ 0 (ARi) - 2J 1 (AR i )j = 0 


(2-24) 


The mode shape can now he -developed with the boundary’ conditions 
expressed in terms of the Bessel’s functions. 

Substituting equation (2-l6) into equation (2-19) for p(r) 
evaluated at r = R^ yields, after simplifying 


A 


Jl’CARi) 


Jl(ARi) 

ARi 


= -B 


Yi’(ARi) 


Y^ARi) 

ARi 


(2-25) 


or 
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A » -B 


[^’ (AR i ) ' — J 

. J^ARi)' 

Jl’(AEi) - — 

L • ARj_ J 


(2-26) 


Thus the mode shape- (2-l6.) eliminating the arbitrary constant A is 

r - - * '■ A ■ 


p(r) = B < Y-^Ar) - 


. ‘ Y 1 (AR i ) 

--55TT 


-Jl’(ARi) - 


JlC-ARi) 


> • 'XRt J 


^l(Ar) > . ('2-27) 


or after use of equations (2-21) and (2-22.) ‘and further simplifying s 

“ * s * » » * 1 

(2-27) becomes • • . 


f • 


P (r) = B ^[ARiJotABi)- - 2Ji(ARi)] Ji(Ar) - (ARiY 0 (ARi) - 2Y 1 (AR i )J J r (Ar)j 


(2-28) 


where 


B = 


B 


ARiJ 0 (ARi)- - 2J 1 (AR i ) 


and since B is arbitrary, we can choose 'it such that B = 1. 
The mode shape normalized at r"= Rq takes the final form 


p(r) = 


|AR i J 0 (AR i ) - 2J 1 (AR i ) 

- 

AR^AR.) - 2Y 1 (AR i ) J-^Ar) 

Jar^A!^) - 2J (AR i )J*Y 1 (AR o ) - 

- 21^) 

K 


(2-29) 

Thus the free-free motion of the elastic propellant segment is 
completely defined by equation ( 2-24) for frequencies and (2-29) for 
mode shapes. 
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Before proceeding to the viscoelastic solution for the propellant, 
brief discussion of three appendices related to the propellant modes is 
in order. General frequency and modal relations, for which equations 
(2-24) and (2-29) are special cases, have been included in Appendix. A. 
Other ‘boundary conditions described by selection of' proper coefficients 
in these general relations are pinned-free and free-pinned shear mode 
conditions. Orthogonality conditions for the free-free modes are 
developed in Appendix B. Also} an alternate method of solving ex- 
pressions (2-24) and (2-29) using asymptotic- approximations to the 
Bessel’s functions is provided in Appendix C. This method of solution 
to the problem at hand was. developed specifically for an early stage 
of the analysis and, to the author 1 s knowledge, does hot exist elsewhere 
in the published literature. It should also be pointed out here that a 

1 f * ; 

computer program has been developed for thi's investigation to solve; 

relations (’2-24) and (2-29) and is included in? Appendix E.< 

* * * * 

» 

Returning to the present .analytical process ^ it is next, required to 

establish the viscoelastic natural motion description, for the pro- 

‘ ' . < . 1 , ’ 

pellant in shear. Equations (2-24) and (2-29) will he used in this 
development. ‘ • 

Viscoelastic Solution 

It is well known that various means can be employed to obtain the 

T. * ♦ 

viscoelastic solution directly from the elastic solution provided 
certain criteria are met such that the dynamical correspondence 
principle applies. Although many quasi-static problems have been 
solved through the use of Laplace transform techniques’, few -classes 
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of dynamical problems have been solved via correspondence principles. 

As is pointed out in the literature [58] , Fourier transform techniques 
have been employed for some dynamical problems,- with* inversion of the 

4 , 

transformed solution after substitution of viscoelastic .moduli often 
being quite difficult* Fisher and Leitman £ 59 ] have developed* a ' 

* * 4 

correspondence principle which utilizes the Fourier transform technique 

t ’ * > ' 4 T 

and which applies to the propellant natural motion being considered 

herein* The authors provide* an analytical proof for the hypothesis 

that if the elastic mode shapes of simple free vibration are dependent 1 

upon only one relevant relaxation function, .the 'viscoelastic mode - 

shapes are the same as those developed for the associated elastic 

*■ * 

material (specifically, those described by equation (2-29)' in the 
present analysis)., and the frequencies are related by the expression 


% 

cue 


2 G(av) 

» 1 + 


G° 


(2-50). 


■where 


cby = viscoelastic natural frequency 

oi e = -elastic natural frequency- established by -equation (2-24) 

G° = the zero frequency, (or elastic) shear modulus 

and 

G(ov) = f G(s) e ~^ s Is 

J 0 

which is the transform of the viscoelastic shear modulus. Equation 
(2-30) evolves from, the following considerations. Using the notation 
of [59], an operator Lq w is> defined as 
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LGS(x,t.) = G°g(x,t) + / G(s)g(x,t-s) ds (2-31) 

J 0 

where G° = the shear modulus at zero frequency and G(s-) is a real- 
valued differentiable function on (0,») whose derivative G(s) is 
absolutely integrable. Also, the integral G(aj v ) is defined as 


G(ov) = f G(s) e " i0 V s ds 
J 0 


( 2 - 32 ) 


For solenoidal motion of a homogeneous isotropic viscoelastic 
solid, the governing equations can be written as 

d 2 u 


-Lq. Vx ( V X u) = 


St 2 


(2-35) 


where u is the displacement function and 


r 

5,' j* , H *. \C * ■ 

G = the shear function 


G° + G(oy) 
•25 


Wow if simple free vibrations are assumed, where the definition of 
simple [ 39) for solenoidal motion is 


V • u- = 0 

then u(x,t) = f(x)e iaiv '^, and by making use of equation ( 2 - 31 ) and 
(2-32), (2—33) becomes 

V 2 f(x) + -^=— - f,(x) = 0 • , (2-3*0 

G° + G(ov) 



16' 


Also for the 618.3X10 regime, if describes a simple free 

vibration, .the resulting equations of motion are 


V 2 ^) + f.(x) = 0 


(2-35) 


Subtracting (2-35) from (2-3*1-) anch simplifying yields . 


CCy 

u“ej 


2 = i + 2bl , 

G° ■ 


. (2-36) 


» T - 

which is the same relation as (2-30) and means that for a simple free 

\ j ’ " • • 

vibration, the mode shapes of 'ihe Viscoelastic solid are the same as 

1 ; . ' 

those of the associated elastic solid 'and the ‘ f requeneies are. related -* 

t * \ * 4 i ' , 

b '/ (2-36). Note that G(<%) is the Fourier transform shear. modulus 
and we seek the value in the frequency domain, such that the ’ 

1 1 

relation (2-36) is usable in its present form. ’ Now it has been 
shown that, for many filled and unfilled' polymers , -the storage and 
loss moduli,. &]_ and G 2 , respectively, may be replaced in -the 
expression G(ccy) = Gq/m^-t- iG 2 (o5 v ) 
by 

Gi = Kq_CD T K (2-37) 

and ‘ • . 


Gp = %n v e 


(2-38) 


where and Kp are constants and k _ and -e -are exponents; 

It has also been shown that for specific frequency ranges of 
interest, the ratio of loss modulus to storage modulus is .quite small,- 
implying small damping. The assumption will be made at this point. 
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therefore, that the loss modulus may he neglected. Therefore, replacing 
G(ov) in expression (2-20) with the storage modulus Gj_ from (2-21), 
one obtains 


i°e. 


2 


Kio>v K 

1 + 

GO 


(2-39) 


Note that frequency dependence of the viscoelastic shear modulus 
is retained with the complex component neglected. This implies that 
the analysis hereafter is viscoelastic in the sense that the shear 
storage modulus is frequency .dependent, hut is quasi-elastic in the 
sense that damping is neglected. Thus, the free-free propellant shear 
motion description is complete, with elastic eigenvalues. A, 
established by relation (2-2t) , and elastic frequencies determined by 



Viscoelastic frequencies, are given by relation (2-39 )j» and viscoelastic 

(as well as elastic) propellant modes are established by equation ( 2 - 29 ) • 

> 

It is now necessary to develop the systems problem by deriving "the 

J r , ►, 

governing equations which couple the Viscoelastic propellant segment 
to the launch vehicle cylinder. This derivation is contained in the > 


following chapter 



III. SYSTEMS ANALYSIS 


To this point, we have developed the propellant free-free normal 
mode solution (2-29) an<l the viscoelastic frequency expression (2-39 ) > 
•both restated here for the purpose of continuity 

- (r) _ C M i J o( M i) ~ - C M i Y o( M i) - J i<^> 

(2-29) 



(2-39> 


These mode and frequency expressions will' be used in the develop- 
ment of the systems analysis which follows. In order to develop the 
governing equations of motion which describe the, system dynamics, we 

will require four basic elements which are as follows: 

* 

(1) Propellant free-free natural modes in shear v 

* x » 

(2) Viscoelastic ‘propellant' frequencies 

(3) Vehicle elastic free-free natural modes and frequencies in 

^ \ 

torsion ' * " . - ' . 

(4) Appropriate physical constraint conditions which couple the 
propellant and vehicle modes together. 

* T ' S 

A graphic description of the system physical components is provided in figure 2. 

The approach used in developing the coupled system problem is the 
convenient Lagrange's equations method, with the equations expressed 
in the -general form for dissipative systems 
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Figure 2.- Launch, vehicle system physical components schematic 
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d/dL 

dt i ^ 


Sl 


+ ^v 


^k. s *k 


\ a l k 


(5-1) 


where L, = the Lagranglan = T (the total system' kinetic energy) — 
•U (the total system potential energy) 

= dissipation function (velocity dependent) 




= Lagrange multipliers 


a^ = coefficients of the Lagrange multipliers 
Dissipation has been considered in order to make the solution more 
general, even though damping has "been considered as negligible in the 
previous propellant dyn ami cs development. 

Now the Lagrangian requires development of the system kinetic and 
potential energies. The kLnetIc energy per unit volume for the 

i 

propellant segment in shear is - 

1 n 

(3-2) 


1 *2 _ 
dT_. = — p d 
*F 2 H rn- 


Now p in (3-2) can be expressed in modal series .form as 

p = ) a” p Ir) 

\ K / . n ^n -l 1 • 

n • 

where 


(3-3) 


. ,• — lot 

• • a =. a e ,• 
n n 


Therefore, integrating over the propellant volume, IrdrdQ, the 
propellant kinetic energy becomes 



■5 18 /' ° r I a n p n (r > 


o ,J R. ■ ^ 
1 . n 


rR -fr~) . ~ 2 

r F = * lS /_ ° r I% P B Cr V 

i \n *' 


Also* the kinetic energy of the launch vehicle can he expressed: -as 

1 f L • P 

• T LV - 2 / z ( x ) e LV U,tr ax * - • 1 ( 3 -! 

. ° 

Again, modal series representation of the angular deflection 0 
in 1 (3-5) will "be used, that is 


where 


9 LV - S °i 5 1 <S) 


— lost 

c. = c. e 
l x 


Therefore., making substitution of (3-6) into (3-5) yields 


T LV = 5 / Z « X C i ? l (s) 


The total system energy is obtained by adding (3-7) to (3-*0> or 


(r) dr + - 


(3-8) 
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Now the potential energy in shear of the propellant per unit volume 
can "be written as 


U. 


F 1- ' 

= o t r. 


vol. 2 rq> rep- 


.. (3-9) 


where the volume is ZrdrdG as defined before, > The shear stress r * * 

; * , rep x 

has been given in relation ( 2 - 7 )" as 

n. ' Sep 

T rcp - G F r 


Also, from relation (2-6), the shear strain 7 ^ is 

Sp 

Substitution of (2-6) and (2-7) into (3-9) an<l integration over the 
volume yields 

a r = i£ 8 *r° G f 


or 


U F ' f'° r f&f " 


(3-10) 


But we have shown in- relation (2-10) that 


-ife- £J\ 

r r 2j 


(3-H) 


and by relation (3-3) 



Therefore 


3F = X a n P n' (r) ( 3 - 12 > 

n 

where, the prime indicates differentiation with respect to r. 

Making substitution of relation (3-3) and (3-12) into (3-H) yields 




(3-13) 


Now substituting (3-13) into (3-10) and simplifying results in 


XL, = itl'G 
F F 


P„(r}' 


/ R ° r I y.' W 1 ‘Z' a n- 
in n 


r 

\ f~ 


dr 


(3-1*0 


which is the final relation for the potential energy. of the propellant. 
The potential energy of the launch vehicle can be written as 

■ v 

* L* d0 T 


U LV = I / J 1T °LT -# (x ' t)2 (3 - B) 


But in (3-6) we had 


LV 


■ X w 


s) 


Therefore 



2k 


^LV 

6x 




B) 


(3-16) 


.Substituting (3~l6) into (3-15} yields 


c 

— 
r— 1 

=2 J LV G LY . 

). h v (s > 

c 

0 

r f 

H-{ 
1 


dx< 


(3-17) 


The total system potential energy is obtained, by 4 adding relations- 
(3-14) and (3-17) yielding . . ' 


U = xlG r 


R. 

1 


L ° r I 1 a n VW -.'I 


n- 


dr 


r L 


L 

z °i £ l' ui 


i 


2 > * * 

dx 


(3-i8> 


Thus the Lagrangian L can he obtained by subtracting (3-18) from 
(3-8) or 

R f L r ~|2 

L = ^ X ° r Z % *n (r) * + \ J Z(x) ' .Z °Z 


i L? 

'R r 


n R -’T— ii r— t P 

J v ° 1 L ** P n' (r) - L\~: 




i L n 

•I* 


- I g lv X, Z “i t 1 * c “). 


dr 


dx 


(3-19) 
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Wow the dissipation function at propellant radius r, assumed to 
"be velocity dependent is . 


c£(r) = i K p(r,t). 2 


■I V 


where K = > k^, the damping coefficient 

■n 

and from (5-5) 


P = I "n 5 n (r > 


n 


or 


p - I a n 


n 


(5-20) 

( 5 - 21 ) 


( 5 - 22 ) 


Substituting (5-21) and (5-22) into (5-20) and integrating over 
the radius of the propellant, the dissipation function c£ becomes 

3= \ L ° X ®ri (r) ] dr (5 ' 25) 


i n 


Now the constraint relations which must be satisfied in coupling 
the propellant to the launch- vehicle at the common interface at R q 
are that displacements and stresses be compatible. The mathematical 
statement of the displacement condition expressed in Lagrange 
multiplier form is 

*1 O’ - B o e Lv] - 0 

Now making substitution of the modal form relations (5-5) an - < 3- (5-6) for 
the propellant and launch vehicle deflections , respectively , yields 



26 


I a n < E o> - E o Z C i 5 1 <S) 


n 


= 0 


( 3 - 25 ) 


Also, the continuity of' stresses across the interface can he expressed 
as 


t r cp “ t r e 

• o LV 


= 0 


(3-26) 


where t it , is obtained from relation (2-7), evaluated at r = R . or 

R CD \ ♦ O 

CT 


rep 




r=R 


r=R 


or 


T R ep = °F R o 1? 


r=R r 


( 3 - 27 ) 


But since modal form series have been used, we can evaluate 
by using, ( 3 - 13 )? or 


3F 


r=R 


Sep 

Br 


R 


r=R 


P„(Rq) 


Z^-Z^ 


n 


n 


(3-28) 


Substituting (3-28). into (3-27). and simplifying yields 


V 


= G. 


F 


z 


a P '(Eo) 
n *n ' 


-l 


P n (Ro), 


.n R 


n 


■o 


(.3-29) 
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Now the shear stress t on the launch vehicle side of the inter- 

A o r tr 

o LV 

face is in equilibrium with the inertia of the launch vehicle, or 


r - E h 05 p 0 T , r 

R o 9 LV s LV 


(3-30) 


By substitution of the modal expression (3-6) for 6^, we 


obtain 


Vnr ' “ 2 ° s X tlU) 


(5-31) 


Substitution of (3-29) and (3-31) into (3-26) yields ” • 

M & F £ a n P n ’(E 0 ) - £ i“ - E o h 05 2 p s , ^(s)j> = .t> 

L n n ' 0 # i J ’ * 

(3-32)' 

T , . ‘ > 

Equations (3-19), (3-23) i (3-25), and (3"32). now provide all the • 

essential elements required to develop the system .governing equations; 

by performing the, indicated operations of (3-l) on these four 

relations in the generalized coordinates q^ = a m , , the * 

following equations are. obtained. For 
>R 


2*15 


r K o fv — 

4 r L "n P n 

i L n 


(r) Pjr) 


m 


- 2rtZG r 


r h 

I 0 


J R. 


/ E ° r X a »V« - I a n- 

• 1 L n - n 

/ k a~ p (r)- 
Z_i n n *n 


pjr) 


n 


■ PJW 


m 


p m (r) to = ^ p m ( Rc ) + ^G p| 


dr 


^ pjR ) 


(3-33) 
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and for q, = c , 


Cj(x) dx 


/ 0 Z<X) I °i 5 i <X> 

- g lv./ 'tvf ^ c i s± M 


^j(x) dx 


- " ¥o^) 


- -h 2 R 0 hof p s 5 ( S ) 


(3-34) 


Relations (3-25), (3-32)', (3-33) and (3-34) completely define the 
coupled systems problem. In order to simplify and render solvable this 
system, recourse to orthogonal tty conditions in the component modes is 
necessary. These orthogonality conditions are derived for the 
propellant cylinder in Appendix B,. and the elastic vehicle cylinder 
orthogonality relations are- also given therein. 

Equation (3-33) can be initially simplified by expanding the 
second integral to 





- r) p m ,ar - / r (v - ¥) 


m 


dr 

(3-35.) 


(where it is understood that the p's are functions of r only and’ 
that the upper and lower integral limits are R Q and 
respectively)'.. Integrating the first of these expanded integrals by 
parts by letting 
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then we have 


R 


{v - rU. i - / il" v to . ; '3-3s> 


V 


* 

The first term vanishes identically at the limits by virtue of the 
free-free boundary conditions (2-17)* 

Now using the D.E. of motion (eq. (2-li ) ) , we have after 
rearranging 


^n 


It 



(- 


0) 


2 5 


n G, 


F 




( 3 - 37 ) 


where the are the viscoelastic frequencies of the propellant 

obtained from ( 2-59) . 

Making this substitution into equation (3-36), we have 



p dr s ^ I— ' / rp P' dr 
m n G_ J ^n 


f p p * 
J m r n 



( 3 - 38 ) 

Now, substituting the right side of relation (3“38) into equation 
(3-35) for the first integral and expanding the second integral, we 


have 
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atl0 F. X a n [“„ 2 ^ />*, n + /l> 1 » V ‘/^ 31 

/ rp p 

p n ' p m ^ + j -t - 2 ‘H (5-39) 

which reduces to 

2irlS co n 2 7^/rp^ dr 
n 

Appendix B relation (B-9) indicates that a value exists 1 only when 
m = n. Therefore,, for all terms of the series when i ^ n, 

21,16 “n 2 I a nJWn * = 0 

n 

and for m = n, 

2*Z5 a) 2 > a / rp p dr =■ a * (3-46). 

n /_. n L / m n / n n T 
“ u n 

n ? , 

, , » ( 

4 I 

where co is the viscoelastic component frequency . > • 

t • f • i 

Also, it can be shown from more orthodox orthogonality , proof s 

* ' * 

i ■ , 

v i 

that the second integral of equation (3~34) reduces to 

" \ 

\ * i : 

> 1 f v 

for i = j ' (3- 

s * . 

« \ f 

where cu^ is the elastic component frequency. 

Now simple harmonic motion in the generalized coordinates and 

has been assumed. Making use of this fact and relations (3-*t-0) .and 
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(3-4l ) , it can "be shown that relations (3053) and (3-3*0 , after 
some simplification, reduce to 


r~~ r> H ” *1 

IT (“» 2 - “O \ + i<4! o 2 X ° k n + >2 S o G F P, 

o * n I 


...,^0)1 

pn ( R o) 


+ R o A 1 P n (R o } = 0 


(3-42) 


and 


'i V K 2 + *) + \ R o M s); + \ R o h “ 2 p s M s) = 0 


(3-43) 

where 3^, = R 2 

n n 

and o> Is the coupled system frequency. 

Recasting equations' (3-42)- and (3-43) to separate system frequency 
elements from natural component frequency elements, 

^R 


h 


ico R 2 f ° k P 2 (r) dr 
oj » n 


a a _ 

ir = \ir% 

o no* 


+ B G. 
o F 


■P n ' (Bo) - 


SlPo) 


B 


O -J 


E o P n ( R o)^l 


(3-tt), 


and 
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r 


3 T. C i + V P s V S) A 2 


L_ 




I V‘ i c i m i " R o \ 


(3-45) 


Also-, the constraint relations (3-25) and. (5-32) can he rearranged 

t icot — icot \ 

(remembering that a n “ a n e and = c^ e ) as follows 


Y, IT R o - E o I °i 


n 


\ " 0 


(3-46) 


and 


Rh p 
o 's 


I °i ^i (s) 


“ - ®p 


IKw-IKw 

* — n u— ' n 


n 


• ' (3-4-7) 

Neglecting the second term of equation (3-Wl*) which is the damping of the 

■* , i , 

nth propellant mode, equations (3-44) through (3-47) can he .expressed \ 

conveniently in matrix form as s 

. * , * 

[M] {Y) =■ M (V), . , . . '('3-48)' 1 - 


0) 


where [M] and [K] are partitioned matrices as .follows . i 





VJ4 


( 3 -^ 9 ) 



and 


H 


*0 


Vo ( V 


E o G F| 








R G T'| 
o pi 


p l (R o } 


*l*<V - 4 


I_ (0 ' 

X 11 


E o p n < E o> 




4 . 


p(R) 

p • R ) - -V 0 .. 

o' R 


V° 




W s > .-Vi (8 > 


0 

1 1 

1 R.p (R ) 1 
. 0 0 0 , 

i.J 

O Fi 

Q t/ E \ p o^ 

ly ( v - — j 

0 

l 1 

Iwv 1 

Vf 

r p >( B j .vv] 


44 tj 

6 

lv>0> I 

J"4v s) j 

i E o°F 

r « pJr H 

p f (R ) - - 4 . -° 
y o' R c J 

0 


I-w b > | 
1 : 


0 

\ s i 2 

• i 

I - E oM s > J 


0 4 

^X(b) I 0 

® * i . 

" 0 

* 

i 1 



0 - ’ 

1 0 1 


0 


VX 

4 =" 


(5-50) 
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and 


R 


R 


n 


R 


( x ) = \ °° 




°± 


^ y 


(3-51)' 


It should, be noted that equation (3-48) represents a coupled 
viscoelastic -elastic- system, because as was stated in equation (-2-39 ‘) 
the propellant frequencies, a> n , are nov viscoelastic frequencies. 
Equation (3-46) is quite general, however, and for the matrices as 
defined, could be either elastic or viscoelastic frequencies. 

It is important to note here one specific advantage of having 
selected the Lagrange Multiplier technique. The use of Lagrange 
Multiplier constraint relations permits the necessary boundary 
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conditions linking the viscoelastic cylinder to "be satisfied on the 1 
stun of the modes, (rather than on individual modes such as is required 
in the Rayleigh-Ritz method) . This artifice allows the use of free- 
free modes and the attendant simple boundary conditions to describe the 
viscoelastic cylinder shear motion. Various slopes of this cylinder at 
the bonding surface can thus be examined parametrically. The effect is 
analagous to permitting a rotational spring to exist at the Interface. 

It should be observed that the order of degeneracy of the systems 
problem matrix is increased by the inclusion of the multipliers, but 
this affords no difficulty to the problem, solution by modern 
computational equipment. 

Solution of System Equations ■ - • 

t , t 

It can be noted that presence of the rigid body, modes and the^ 

T / , ‘ 

Lagrange multiplier relation causes singularity of both matrices [N] 

? T 

and [M]', respectively, thus precluding solution of the eigenvalue-, 
problem in its present form. ' An. alternate, form [-reference -4-1} for the c 

« f i - 

problem is obtained by scaling the eigenvalues as 'follows: 

Equation .(3-48) can be written as 

[M - M + [t|»] - [tjMjJ {y} = 0 ’ ‘ ' (3-52) 

where y is a constant scaling factor. 

Factoring and rearranging relation (3-52) yields 

[M] (V> = -- 1 -p I* + T]M] (y) = 0 

^ tj + or ^ J 


(3-53) 
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which is valid for 


> -0) 


The square matrix on the right-hand side of equation (3-53) is now 
nonsingular, and equation (3-53) can he- written as 


[U] (y) = V {y> 


(3-5*0 


where 


[H] = [N + nM]" 1 [M] and V = — 1 g - 

• T1 + CD 

Equation (3-5*0 is now in standard eigenvalue form and is solvable by 
general iterative procedures and sweeping methods (outlined in 

1 ‘ 

Appendix D) for the eigenvectors and respective natural frequencies-. 

Obvious advantages to this method are that the solution yields 
all inodes, including the rigid body modes, without recourse' to sweeping 
matrices for establishing the lowest system fundamental. There are ’ 

t * * 

minor difficulties to be noted in this method of scaling eigenvalues, 

t 

however. First, lacking prior knowledge as to the lowest frequency, 

-* ^ i 

q must be chosen as positive. With positive,- q .close proximity of - 
adjacent latent roots forces slower convergence of the solution than 
would otherwise occur in the unmodified system. Also, the order -of the 
[H]' matrix is larger than required since the rigid body mode's, are. 
present in the modified form. These minor disadvantages are over- 
shadowed, however, by the fact that the form of equation (3-5*0 is 
quite familiar and is easily solved by elementary means . 
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It is readily recognized that the use of digital computer solution • 
•methods are required for expedient solution of matrix expressions such 
as ('3-5*0 ♦ A computer program was developed specifically for the 
problem described herein and has been supplied in Appendix E.' 

It will be noted that the eigenvector as defined by relation 
(3-5l) represents the modal coefficients of the ' propellant and vehicle 
natural modes (also included in the eigenvector are the', constraint . 
forces, Therefore in order to obtain the deflection characteristics 

t 

of the components at any system frequency, it is necessary to multiply 
the resulting eigenvalue, <^Xy>, by a matrix comprised of the component 
natural modes in the maimer following: By selecting the partitioned 

eigenvector section relevant to propellant modes, that is 



and using the natural modes of the propellant, p n (r) , described by 
relation (£-29) with the radial coordinate given by 
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( 5 - 55 ) 


where ^VpJ> represents the propellant .deflection at the qth 
system frequency. A similar matrix operation which will yield the 
launch vehicle system deflections is 



4o 



( 3 - 56 ) 


where the longitudinal coordinate is described, by 

j — 1. 2 * * * » k 

J 

and ^y represents the vehicle deflection at the qth- system 

frequency. The £ j(x) are not developed in the present investigation 

since methods for establishing natural mode data, for launch vehicles 

in torsion are veil documented ['^2] .- 

The total systems problem has now been cast into matrix form, and 

. - . ' 

appropriate solution methods have been developed. The - next logical 
step in the investigation was to select a satisfactory means of total 
•systems analysis verification. This process is described in the 


following chapter. 



IV. VERIFICATION OF SYSTEMS ANALYSIS 


As lias .teen stated, the primary goal in the investigation was to de- 
velop a systems analysis which could be applied to a characteristic 
launch vehicle such as that for Project FIRE shown schematically in 
figure 2. Prior to application of the analysis to the FIRE system, how- 
ever, it was deemed important .to use an experimental model to verify the 
fact that the analytical method would correctly predict response . 
characteristics of a system comprised of two components when the modal 
information for the components had been analytically developed separately 
•and then the components coupled together. 

While it would have been very desirable to employ a large-scale 
dynamic model of a launch vehicle with a solid propellant segment such 
as that discussed in Chapter I, an appraisal of fabrication and test 
setup costs required to develop such a model indicated the unfeasibility 
of doing so. After giving careful consideration, to several simpler 
vehicle system models which could be manufactured, and after noting that 
the frequency response range of the cylindrical vehicle representations 
would be unacceptably high, it was obvious that other avenues of experi- 
mental verification would have to be explored. It was then observed • 

' j - * 

that in the analytical development of the two boundary conditions for the 

I 

system (equations ( 3 - 25 ) and ( 3 - 32 ) ) , the radius of attachment could be 

h 

changed from the outside to the inside radius with little effort. ‘ - 

.1 * 

Although the resulting experimental model then bore little similarity to 

„ . < * < 

a launch vehicle system, it was concluded that a neoprene disc could he 

* ♦ I f 

attached at its inner radius to a brass shaft and by proper selection, of 
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physical characteristics such as radii and lengths., composite system 
frequencies would- he in a sufficiently low range (below 1000 cps) as to 

he readily measurable. Such a model was subsequently constructed and is 

/ « 

described schematically in figure 3. In addition to the neoprene disc 
fabrication, a steel disc, equal in rigid body inertia to the neoprene 

* i 

item, was made. Use of the steel disc will be discussed subsequently. 
Physical data for all system components are given in Table I. 

■ 4 , i 

Disc modes and frequencies for both the neoprene and steel items 
were established using the propellant shear mode computer program' of 
Appendix E which was developed to solve equations (-2-24) and ('2-29) of 
Chapter IX. Elastic assumptions for the neoprene disc were used 'because 
it was expected to behave elastically for small-time, short-period 
frequencies. Shaft modal data were determined by a standard torsional 
modes analysis not included herein. Coupled neoprene disn-brass shaft 
and steel - disc-brass shaft frequencies and 4 modes were then established 
via the systems analysis computer program of Appendix: E which was written 
to- solve the matrix relation (,3-54). The ana3ytical free-free mode 
shapes of the neoprene disc have been presented in figure 4, and the 
resulting analytical neoprene disc deflections at system frequencies are 
shown in figure 5. Steel disc analytical modal data were not plotted, 
since the more salient information for the steel disc-shaft combination 
was the frequency spectrum which is discussed later. 

Once the analytical information had been established, an experi- 
mental program was undertaken to verify the adequacy of the systems 
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Figure 5»- Disc -shaft experimental model schematic 


TABLE I. - SHAFT-DISC PHYSICAL PARAMETERS 


Item 

Material 

Dimensions 

Remarks 

Shaft 

Red brass 

1-i'n. nominal diameter 
100- in. length 

Solid rod with two 
l/2-in. -thick by 
2-in. -diameter 
brass collars for 
disc attachment 
(see sketch) 

"Flex" 

disc 

Neoprene 

Ik-in. outside diameter 
1-in. inside diameter 
1-in. thick 

50 durometer 

Shear modulus = 
90 lb/in^ 
nominal 

"Rigid" 

disc 

Steel 

8.6-in. outside diameter- 
1-in. inside diameter 
1-in. thick 

- 



Brass shaft configuration 






















Figure 5 «- Analytical neoprene disc characteristics at system frequencies 


analysis. The test setup and apparatus are described in the next 
section. 


Test Method 

The boundary conditions simulated for the system were free-free. 

A small cable sized to ensure negligible torsional restraint was used 
as the upper end suspension method (-fig. 6). The lower end was com- 
pletely free. Two ’five-pound shakers (fig. 7) were hooked up l80°' 
out of phase to provide as clean a torsional input as possible to the 
shaft. The driving point was near a system node' such that very small 
forcing was required to excite system modes. To ensure that only system 
torsional response characteristics would be measured,, a complete trans- 
verse modal survey of the shaft was' performed prior to the torsional 
test to identify other modes. Also, several rotational orientations of 
the accelerometer mount disc radius were evaluated to assess the possi- ' 
bility of spurious mode detection. Minor changes in response measure- 
ment were noted after changing the disc position, and the tests were 
completed with the accelerometers aligned perpendicular to the back- 
board (fig-. 8). 

The neoprene disc (fig. 9) was supported prior to tasting by a ply- 
wood disc covered with teflon to preclude neoprene disc sag (fig. 10). 
The disc was fixed between two brass collars and double-backed tape 
between collar and disc on both sides- to preclude relative rotation at. 


the interface. 



Figure 6.- Free-free suspension for experiment. 
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ACCELEROMETER ORIENTATION I 
TOP VIEW 



ACCELEROMETER ORIENTATION II 
TOP VIEW 


Figure 8.- Radial accelerometer orientation schematic 
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Figure 10.- Neoprene disc location in system 
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Measurement data were obtained from five small (5g) Unholtz -Dickie 
accelerometers mounted across a radius (see fig. 10). All accelerometer 
readings at each resonant frequency were recorded on an oscillograph. 

Five runs were made at each resonance allowing resonance verification 
via oscilloscope monitoring of a separate accelerometer Lissajous 
pattern for each run. Frequencies of the resonances were read from a 
digital counter and from the control oscillator. Final determination of 
frequencies was made later by direct reading of the oscillograph traces. 

The complete test setup for the neoprene disc-shaft is shown in 
figure 11- The setup was identical for subsequent tests with the neoprene 
disc replaced with the steel disc (fig. 12). For the steel disc modes, 
only three accelerometers were used. 

Test Results 

All data are presented in Table II for both shaft-disc systems. The 
analytical data generated via the systems analysis method developed in 
Chapter III and discussed in the prior section are also presented for 
comparison. 

The experimental and analytical neoprene disc deflections at the 
disc-shaft first mode frequency are plotted in figure 13 . Frequency 
agreement was generally acceptable for the experimental setup utilized, 
but corroboration of the analytical mode shapes was difficult. Several 
considerations concerning the experimental mode shapes are as follows: 

1. Only five accelerometers were used, making node determination 
at higher frequencies (where nodes are more closely spaced) difficult. 



Figure 11.- Complete test setup 
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TABLE II. - RESULTS OF SHAFT-DISC EXPERIMENT AND ANALYSIS 


System 

•Mode . 

Frequency , cps . ' 

1 « i * , 

1 Analytical 


Shaft - 
neoprene 
disc 

1 


'10? ‘ ‘ 
- f 

2 

162 

* • 162 

3 

1 

250 . 

k 

286 

27.8 

5 

3bQ 

350 

6 

399 

1^24 

Shaft - 
steel 
disc 

1 

266. 5 

266 

mm 

| 

* 


*Not included in test equipment frequency span 
range. 
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2. Close proximity of several lateral and torsional natural 
frequencies also hampered mode determination at higher frequencies . 

3- The efforts to obtain good fixity .at the disc-shaft interface 
may have affected the modal responses adversely in that there is the 
possibility that the tightening of the collars may have squeezed the 
neoprene disc at the inside radius thereby shifting the effective 
attachment radius outward. Several analytical runs with different 
attachment radii were made by computer in an attempt to achieve better 
•modal agreement in the higher modes, but better results were not obtained. 

With these shortcomings, however, it is felt that the systems, analysis 
predicted the frequency responses well. It should he noted that the 
frequency spectrum for the rigid steel disc (even though possessing rigid 
body inertia identical to the neoprene item) differed considerably from 
that for the neoprene disc thus pointing up the errors to be encountered 
when flexible body effects are neglected. 

Having established what was felt to he satisfactory agreement between 
analytically predicted and experimentally determined natural system 
characteristics, particularly in the frequency spectrum, the final step- 
in the investigation was to apply the analytical method to the FIRE 
vehicle. Details of the example problem are given in the following 
chapter . 



V. EXAMPLE SYSTEMS .PROBLEM 
‘ # ) * 

v 1 * 

The systems analysis developed in Chapter HI is applicable to the 1 

' ' ■ ■ ' • v ; * 

examination of a wide array of composite systems torsional dynamics 
problems, including systems such as tbe experimental model discussed 
previously. However , it was developed for the .specific purpose of 
describing the torsional -motion of a charact eristic launch vehicle such 
as used for project FIRE. The vehicle system-, is- comprised of a near 
axisymmetric discontinuous thin-walled elastic cylinder with a solid 
propellant motor attached near the forward end in the payload region 
(fi-g. l4). The case examined is not for academic purposes only, it may 
be noted, • since flight data have 'indicated that- the three booster 
engines which are gimbaled to permit roll control as well -as thrust, 
tend to "chug" on shutdown in a. manner which induces a quasi -sinusoidal 
torsional excitation (fig-. 15)- devoid of any components- which" would 
induce other modes of vibration. This input is in turn transmitted 
forward through the length of the vehicle, ultimately reaching the 
propellant attachment station. Since the Young's modulus for a typical 
propellant casing is E = 0.6 X 10^, a nominal zero frequency value for 
the propellant is E = 500-, and the propellant -length fraction of the 
total system length is small (approximately 6 percent), classical 
torsional deflection over the casing length will be neglected for the 
application of the analytical method. It will also be noted that 'for 
the rigid casing assumption, the case responds by undergoing angular 
rotation constant throughout its length under the influence of 
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torsional excitation. Also, when the difference in casing and 
propellant moduli as well as the mounting arrangement are considered, 
justification for the assumption that the propellant would "behave in 
pure shear i-s evident. Figure 16 contains more detailed information 
about the payload portion and is typical of the mounting arrangement 
for flight research experiments of this class. The viscoelastic 
propellant is shown shaded and is encased in a- fiberglass -wound casing 
(item l) . The load path in twist,. is through a bolted joint near the 
propellant aft end, via the aluminum forging (item 2) through the "V" 
clamp joints in the adapter section, then through the aluminum skin 
structure to the mating joint (item 5)- The sandwich structure (item 3) 
is not intended as a load path and is for the purpose of flight 
accelerometer mounting only. 

Discrete Vehicle Model 

The mathematical model used to describe the longitudinally 

discontinuous launch vehicle is given in figure 17 .. The physical 

attachment station of the propellant is as shown in the figure. The 

discrete torsional stiffness and inertias are given- in Table III. It 

may he noted that for the vehicle modal deflections and associated 

generalized inertias computed by standard eigenvalue routines using 

* 

the discrete data, only the modal deflections of the attachment station 

( ' t 

(i.e., the launch vehicle -propellant mechanical interface) are required' 
in the systems problem. ('Also, the propellant container inertia is 
lumped at the vehicle attachment point and considered in 'the* vehicle 
modal data . ) 



MATERIALS 


FIBERGLASS 

2014-5652 FORGING 

7075-T6 A1 - MAHOGANY SANDWICH 

2024-T3 A1 SKIN 

2014-T652 FORGING 


/ 


"V n CLAMP JOINT 
BOLTED JOINT 


Figure l6.- Splid propellant mounting schematic 
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•ROTATIONAL INERTIAS 
ROTATIONAL STIFFNESSES 



(DISCRETE BRANCH MODEL) 


Figure 17.- Atlas math model, 












TABLE III.- TORSIONAL PHYSICAL DATA OP 




f 

I Roll 

X lO*^ 



wm 

1 lb-in. ^ 

“PTm V 1 0*^9 




- 

iup A _LVJ 

Node 


X-259 with 

x-259 With- 

in. -lb/rad 


mu 

solid fuel 

out fuel 

■19 

79Q 

_ 


_ ^ 

25 

26 

8o8 1 
830 

5-215 

5.215 

5.215 

5.215 

8.660 

3-.2.30 

27 

842 

9.310 

9.310 

7.640 

28 

872 

32.4oo 

32.400 

4.290 

29 

889 

8.230 

1.84o 

5.224 

50 

894 

7.600 

7.600 

.17.762 

31 

908 

-9.730 

9.730 

2.214 

32 

926 

.4.100 

4.100 

2*278 

33 

954 

5.910 

5.910 

1.045 

34 

996 

.860 

.860 

.213 

35 

976 

1.925 

1.925 

— 

36 

948 

6.24o 

1.840 

.0681 

37 

930 

10.220 

2/660 

.167 

58 

910 

( 10.220 . • 

2.660 

.150 

39 

870 

.671 ‘ 

. .671 


40 

842 

.459 ' 

.459 ’ 

.0143 

, 29 

889 

— 

— * , 

— — ^ > 

38 

910 

-- 

-- 

■ - .143 

29 

889 


— 


39 

870 

— 

! 

v .OI83. 
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The launch vehicle torsional modes and frequencies with the 
propellant excluded are given in figure 18 and differ slightly from 
those of reference [42], in that the propellant casing is considered 
rigid in this study. 

Propellant Model 

Recourse to discrete analysis is not required for the propellant 
since it can he represented hy the continuous deflection result in 
shear (developed hy the closed form solution in Chapter II for a visco- 
elastic thick walled cylinder. Some license was taken, however, in 
permitting the mathematical model to 'have a circular cross-section 
port. In reality this inside port deviates slightly from circularity 
although it has heen assumed that for this study the system 

dynamic differences due to the assumed and real cross-sections are • 

■ 

* 4 

negligible. J ' * 

Ehysical data for the propellant » are provided in Table XV > 

- * 

Vehicle Systems Analysis * - * , ** 

i i * * * 

The method employed in solving the systems problem was to follow , 

i * -t ' 

the same analytical steps as described in Chapters II and III. In 

- * * \ 

order to evaluate various propellant assumptions, shear mode' shapes 

with frequencies for an elastic propellant (fig. 19 ).were first • . , 

established using equations ( 2 - 29 ) and -(2-24) and then combined with 

* » . * * 

launch vehicle torsional modes to develop the systems frequencies. In 
the systems analysis for this propellant assumption, 9 propellant modes 
and 12 vehicle modes were used. Modal and frequency results of the 
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Figure 18. - Concluded 
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TABLE IV. - PROPELLANT PHYSICAL DATA 

Outside radius 14.85 in. 

Inside radius 4.5 in. nominal 

Length 5 6 in. nominal 

Density .0635 lb /in? 

Shear modulus 200 lb/in 2 nominal 














1362 rad/sec 
MODE 4 



406 rad/sec 
MODE 1 




Figure 19 i- Propellant shear mpde shapes 
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1 

coupled system dynamics analysis with the elastic propellant are 

* " J 

presented in figure 20 for the propellant and figure 21 for the 

“ ‘ t 

vehicle . , , • : 


Mext the viscoelastic propellant was considered with frequencies 

i ' , 

computed from equation (2-39) , restated here' for clarity 



(• 2 - 39 ) 


Constants were obtained from modulus versus frequency data of 
reference (21 over the frequency range of 300 bo 3000 cps. It was 
clear that for the storage modulus as herein defined, that is, 



K.0) K 
1 v 


selection of = 2.^7 and k = 1 would appropriately match the 
referenced data for the frequency range mentioned above . 

o / 2 

Propellant manufacturer’s data. yielded a value of G - 200 lb /in 
(G° = the shear modulus at zero frequency). 

Using the above data, viscoelastic frequencies were computed for 
the first five propellant modes and the propellant characteristics 
again coupled to the vehicle modes through the systems equations. The 
resulting propellant deflections for the -viscoelastic assumption and 
the associated vehicle modes at the system frequencies are presented 
in figures 22 and 23 , respectively. 



0 ) in rad/sec 
Ro - 14.85" 
Ri = 4.5" 



Figure 20.- Elastic propellant shear deflections at system frequencies. 
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to = 297 rad/sec 
= 322 rad/sec 
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Figure 21.- Launch vehicle deflections at system frequencies with 
elastic propellant. 
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Figure 23.- Launch vehicle deflections at system frequencies with 
viscoelastic propellant. 
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- ' < * * 

Finally, the propellant was assumed rigid ( and the systems analysis 

! 

again performed to yield the vehicle deflections at the system 
frequencies (fig. 24). Because of the simplicity of the propellant 
assumption, no propellant modes are presented. The results agree 
well with those of reference [42] even though the method of analysis is 
quite different . 

For comparison, the system frequency spectrum for each of the 
three propellant assumptions ' is- presented in- Table V. 

Discussion of Results 

It is readily apparent from observing the Table V- results that 
there are essentially no differences in the frequency spectrum between 
the rigid and viscoelastic propellant -vehicle systems. The differences 
between the elastic and viscoelastic frequency spectra are more 
pronounced with- six. elastic propellant-system modes occurring in the 
same range (0 to 696 rad/sec) as -that in which only four viscoelastic 
propellant -system modes occur. The mode-by-mode frequency equivalence 
in the viscoelastic versus rigid propellant vehicle systems results, 
however, are misleading. Whereas the viscoelastic propellant behaves 
rigidly through the third system mode (511 rad/sec) (see fig. 20)' it 
will be noted that the propellant responds at 696 rad/sec in quite a' 
flexible manner, 

A further observation may be made by comparing the system vehicle 
response data of figures 21 and 22. It is noted that very little 
differences occur in these modes with viscoelastic versus rigid 
propellant assumptions. It appears that the rigid propellant assumption 





VEHICLE 

TORSIONAL 


i - 



0 200 400 600 800 1000 


x. Inches 

Figure 2h-. - Launch vehicle deflections at system frequencies with 
rigid propellant. 
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TABLE V.- LAUNCH VEHICLE SYSTEM MODES AND FREQUENCIES 
FOR VARIOUS PROPELLANT ASSUMPTIONS 


Mode 


2 


3 
■ b 


5 

6 


•Launch vehicle system frequencies (: 

rad/sec) 

Elastic propellant 

Viscoelastic propellant 

Rigid propellant 

297 

252 

252 

322 

382 

385 

4 76 

511 

511 

502 

696 

'696 

588 

7^2 

7^2 

696 . 

886 

886 
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has small effect on the vehicle responses in a gross sense as compared 
to the viscoelastic propellant assumption. The appropriate propellant 
assumption however, is still quite important since it is always 
necessary to have complete knowledge of key system component natural 
response characteristics to insure structural design integrity as well 
as proper guidance system selection. 



VI. CONCLUDING EEMAEKS 




The- primary objective of the enclosed' analysis was to develop a 

t 

systems analytical method which would permit the examination of various 

\ „ 

upper stage propellant assumptions and the effect of these assumptions 
on launch vehicle system dynamics in torsion. As a result of the 
investigation and in order to meet the objective, a launch vehicle 
systems analysis has been performed wherein the propellant and vehicle 
natural modes were generated separately and then coupled together via 
Lagrange multiplier relations for the appropriate boundary conditions. 
For the propellant shear displacement., a general closed form Bessel’s 
solution to the governing equations of motion has been developed. 

Simple harmonic motion was assumed, and general expressions for both 
the secular equation and the propellant normalized modes have been 
derived. Also, an- orthogonality proof for the propellant free-free 
modes has been provided. An alternate and concise method for 
evaluating the normalized free-free propellant modes and frequency 
expressions using asymptotic approximations to the Bessel’s functions 
has-been presented. 

The viscoelastic frequencies of the propellant were developed by 
neglecting loss modulus and using an empirically established frequency- 
dependent storage modulus in conjunction with an existing frequency 
correspondence principle relating viscoelastic and elastic frequencies. 
Since damping was -neglected, the analysis may be considered as quasi- 
elastic but with variable shear modulus . Assumptions that damping was 


8l 
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negligible were justified on the basis of published data on' the 
relation of storage modulus to loss modulus for the frequency range of 
interest ., 

In the systems problem development., several convenient mathemat- 
ical methods were employed which rendered the problem more tractable . 
First, Lagrange multiplier methods were used to describe the constraint 
condition s,. Boundary conditions could thereby be satisfied on the sum 
of the vehicle and propellant modes rather than on individual modes 
such as required in Raylexgh-Ritz procedures. Readily* obtained free- . 
free component modes can therefore be used. , Also., use of Lagrange , 
multipliers causes convergence on the exact frequency from ~bhe lower 
side, as compared to Rayleigh -Ritz procedures which converge frcjm the 

i 

’higher side. Although use of the multiplier relations uncouples ^the 

, * 

* 4 * \ t *■ 

equations of motion for the vehicle and propellant,, the -effect of the 
multipliers is to increase the size of and cause singularity of the 
system matrices. This problem is alleviated by application- of a 

J- 

frequency scaling technique which, permits even rigid, body elements in 
the system matrices. * 

Application of the total analytical process was effected in two 
ways. First, an- experimental model consisting of a neoprene disc 
attached at its inher radius to a brass shaft was built to verify the 
analytical method. The shaft and neoprene disc analytical modes were- 
generated separately, then coupled together in the systems problem. 

The experimental frequencies and mode shapes were then obtained and 

r 

compared with analytical values. Good agreement was obtained in the 
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frequency spectrum although only the lower modes could he obtained 
experimentally. It is felt that’ more accelerometers should have been 
used, or perhaps other methods should be developed for modal 
determination of such systems. 

The systems analysis, was. next applied to a characteristic launch 
vehicle and rigid, elastic, and' viscoelastic propellant assumptions 
were examined. . It was found that -while the vehicle modal responses at 
system frequencies differed only slightly for the three assumptions, 
propellant component natural frequency and responses at system natural 
frequencies were considerably different. It may be concluded that 
correct propellant description is more important from -the standpoint of 
its own response at system frequencies rather than from any modifying 
effect it may have on overall vehicle response. 

If is also concluded, that appropriate propellant flexibility 
assumptions are required if correct response characteristics of launch 
vehicle systems which contain solid propellants as upper stages are 
to he determined. In this regard, it is felt that the dissertation 
investigation extends previous work in vehicle torsional dynamics. 

A more formidable set of tasks which are recommended to extend the 

sT ' 

» * * 
analysis contained herein would include* three-dimensional visco- 
elasticity solutions (which would thereby permit torsional as well as 

r * 

shear deflection throughout the propellant .length), elastic propellant 

‘ _ \ 

case assumptions and the use of noncircular propellant core 


configurations 
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IX. APPENDICES 


In this chapter are included studies relevant to the foregoing 
investigation. 


Appendix A - General Frequency and Modal Relations 
for Propellant Shear Dynamics 


The general solution to the Bessel’s equation (2-15) describing 
the propellant motion in shear (2-16) was found to he 

p(r) = AJ^r) + BY^r) (A-l) 


Free -free boundary conditions were used in Chapter II to establish the 
natural modes for' the investigation. The resulting, equation for the 
frequencies (2-24) was 


IlWcW - 2 j i<* r o 0 CwoO-V - 2 
- C\W»B 0 > - 2 - 2 = ° 

(A-2-) 

and the expression for- the normalized modes- (2-29) was established as 
- 2 J-^XR^Y^Xr) -'Qd^XR.) - 2 J^Xr) 

C^oC^) - 2 J l (AR iD Y i( AR o ) - [^(XR^ - 2 J l (XR o ) 


Pinned-Free Modes 


(A-3) 


For the case of pinned-free boundaries of the propellant, the 


boundary conditions are expressed as 


89 
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p(r) 


= 0, 


R i 


ap 



p(Bb) 

Ro 


(A-4) 


Making substitution of relation (A-l) into (A-4) and performing the 
indicated operations yields 


and 


AJ 1 (AR i ) + BT 1 (.XR ± ) = 0 


AJ r (XR ) + BY (TsH ) 

"i'O'V + - R — - 

o 


(A-5) 


(A - 6 ) 


from which, the transcendental equation is readily developed by 

J , 

rearranging and setting .the determinant of the coefficients equal to 
zero, or 


J 1 (AR i ) 

Jl'(AR 0 ) 


J, (*R J 

1 o 


V (jE o> 




- 0 


; (A-7) 


, - «■ * 

Making use of the recurrence relations for the Bessel's functions 1 
(’2-2-1) and (2-22) permits (A-7) to be simplified to 


(A-8) 

- - 2 .C J i (7B iO - 0 
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the normalized inodes expression is developed in. ,a straightforward 
manner to yield 


5Cr) - 


[?i (XR iO -[ 

Y i (AR iO J i (7a>) 

&.(**! 0 Y l (XR o ) ■ 


] J-lCBo) 


(A-9) 


where the normalizing radius is R o . 


Free -Pinned Modes 

For the case of free -pinned, modes, the boundary conditions of (A-6) 
are readily converted simply by interchanging quantities R q and R^, 
or 


aj^xr^) + by 1 (xr q ) = 0 

AJ (XR ,)• + BY _(XR. ) 

AJ-^XR..). + BY 1 , (XR i ) = — i-g- — (A-10.) 

which leads directly to the transcendental relation 


J lW 

^’(xr^ 


J-^XR^ 


R. 


r i< W o> 


T, (W, ) 

Y i^i> 'V- 


(A-ll) 


The normalized modes relation is readily expressed as 


jAR i J 0 (XR ± ) - 2 J-^XR^] Y-^Xr) - ' 

[AR i Y 0 (XR i ) - -2 Y^XR.)" 

|. Jl(Xr) 

JxR^oC^). - 2 Y 1 (AEo)- 

“xR i Y 0 (XR i ) - 2 Y^(XB^)" 

1 


* (A-12) 
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By examination -of relations (A-2), (A-3), (A-8) , (A-9.), (A-ll) , and 
(A-12) , several observations can be made. 

First the transcendental relation can be expressed in general 
form as 

■ &i ’A J 0 ( »V - a s C a 3 - a 4- Y i <w, iO 

- C a 'l »o V’V - a 2 V^oO C a 3 J 0<*V - a 4 V M iO - 0 

* 

. ' . (A-13) ' 

f i 

and the normalized modes equations can ,be expressed in general form 
as 

_ [a^.J 0 (XRi) V^^j] -.^.Y^^^-a^C^jJ^CAr) 

^AH.JqCXR.) a 4 J L ( ja^Yo.tAR.) - Vl^ilK^ 

: '* ’(A-14). 

*- - f 

where by appropriate selection of the a coefficients the three cases of'* 

* f 4 

free-free, free-pinned, and pinned-free modal characteristics can be * 
established. The proper coefficients are given in the following list. 



, CASES 



Free -Pinned 

Pinned-Free ! 

a l 

1 

0 

1 

a 2 

2 ~ 

-1 

2 

a 3 

1 

1 

0 

% 

2 

2 

-1 
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An additional observation is that the mode shapes (A-l4) for the 
free -free and free -pinned boundary conditions have identical form with 
= 1 and = 2 for both cases. The identical nature of the modes 
is due to the selection of the free end condition for both cases to 

i 

develop the modal expressions. This selection is clarified by the 
steps required to proceed from relations (2-25) to (2-29) in Chapter II. 
The transcendental equations for the two cases are of course different. 

A general set of equations for free -free, free -pinned, and pinned- 
free- modes and frequencies has thus been developed in this investigation, 
with free-free modal information being used in. the earlier systems 
analysis. 



Appendix B - Orthogonality of the Component Modes 


Thick -Walled. Cylinder Component 
The differential equation of motion for the element of cross 
section is (from equation (2-l4)) for any mode 1 


a £ Pn 

dr 2 


1 a Pn ./(«n 2 S 1 \ 

r dr \ Op " r2 ] Pn 7 


(B-l) 


where p n indicates the normalized nth mode p n (r). * . 

■f ' 

The same equation written for the normalized rath mode, p^r) is 


a 2 p m ! dp m . K 2 B 1 \ 

+ - — + I r 1 P m = 0 


dr' 


,2- r d-r 


\ % r 2 / 


(B-2) 


Multiplying (B-l) by p m and (B-2) by p n , subtracting the two 
results and simplifying yields ' * 

d 2 p, 


0 


,d 2 p. 


dp. 


("V 2 - j rltaPn - r Ph— +Pn-te' 


&r c 


-n 

dr^ 




'm „ o ~ 9m 


n 

dr 


or 


Gy 


(con 2 - c%2) rpm p n = p n ± (rp^) - Pm (rp n ’) (B-3) 


Integrating both sides of equation (B-3) Between the limits R^ 
to Rq yields 

1(^2 - -lift ^ -/ b _° P„ (Uta'l to - j^° ’a, i(rp n ') dr 
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The right side of equation (B-4-) can be evaluated by integrating 
the first of the integrals by parts* letting 


u = p n dv = teW) dr 

du = p n * dr v = rp m * 


Thus, the first integral on the right side of equation (B-4) becomes 




miPm' 


■>JEL 




~ f rPm'Vn dr 

J B ± 


(B-5) 


Similarly, the second integral becomes 

*o 


E PmPn' 


■Bo 


R-? 


^VV dr 


(B -6) 


Substitution of equations (B-5) and (B-6) into equation and 

evaluation of the result at the integral limits yields, after 
simplification . . 


—• (con 2 - <% 2 j J ° rpjjjPn dr = R 0 P n (Ro)pm’ (^o) “ RoPm(Ro)Pn’ (Bo) 

J ' E i 

- RiPn(®i ) Pm ' 0% ) + BiPj ft (R^)p^ , (Ri) 


(B-7) 


Substituting into relation (B-7),. the following free -free boundary 


conditions 
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P n '(Eo) 

Pn r (Ri) 


P n ('Ro) 


P n (Ri) 

~^r~ 


which are true for all inodes, it is seen- that the stun of the right-hand 
terms of (B-7) vanishes, or 



and, since the o> n are distinct, then for the total cylinder volume, 

4 ; j 1 

the orthogonality condition is 


‘ 2**5 


*PmP n dr = 0 , m 7 ^ n 1 - 


Ri 


= Mc> , m - n 
^n 


where Mji^ is called the generalized mass of the fuel. 



Elastic Thin-Walled Cylinder Component 
It can readily be shown from established procedures that a like 
orthogonality relation holds for the thin-walled cylinder. The 
expression is 



(B-10)’ 
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Appendix C - Asymptotic Approximations 


Unless Bessel’s function subroutines are readily available, 
computation of the generalized mass (equation B~9) in the orthogonality 
derivation requires cumbersome evaluation of an integral in which the 
Bessel’s functions must be expanded in series form. It is convenient 
to utilize asymptotic approximations to the Bessel’s functions in that 
computation. The computation of the generalized mass and the 
associated reevaluation of equations (2-2^) and (2-29) follows: 

The expression for the effective mass of the fuel as developed in 
Appendix B is 


Mj = 2*26 


/J 0 r Pn( r )j d ar 


(c-i)’ 


where 


p n ( r ) 


P n (r) is obtained from equation (2-29) as 


J^7vR i Jg(AR i ) - 2J 1 (A%) 

Y L (to) - [ 

^i Y 0 ( AR i) - 2Yi(XR i ) J x (Ar) 

^JgCTdRi) - 21^^) 

Yl'OMU - 

|\Rj_Yo( A^i) - 2Y^(7\Eq)^j Jq(ARo) 


Substitution of the right side of (2-29) for the normalized mode 
shape P n (r) into equation (C-l) yields 




r R o jLARiJoCARi) - 2J 1 (AR i )]Yi(Ar) - [ARjYo(ARi) - 2Y 1 (AR i )j J x (Ar) ( 2 
Mp - 2itl5 / r < . . -> - r _ ) 

' *i [^iJ 0 (ARi) - ; 2J 1 (AR i )jY 1 (AH 0 ) - ARiYoCARi) - SYiCARi) J^ARo)! 


dr (C-2) 


where A = A n the eigenvalue of the nth mode. 

* » 

Using the asymptotic approximations for Bessel’s functions of the first and second kinds, 
which are: 


J n (x) = 


„„„ /— it nrt\ . 

cos ( x “ F ' T) 


(0-3) 


Y n (x) a 


• / rr 3QTC i 

"■ i z T/ 



(0-4) 


and substituting equations (C-3) and (C-4) into equation (C-2) and simplifying, the numerator of 
the braced quantity becomes 

' |ARi.sin (Ar - cos (AR t - £) - ARj_ cos (Ar - sin (ARi - J-) ]| 

S'' 


r 


< 


,jtA \JR±t 


AR^ sin ^Ar - cos (l\R± - - ARj_ cos ^Ar - sin ^XR^ - 

[+ 2 sin ^ARi - ^ cos ^Ar - - 2 cos ^ sin ^Ar - 


(C-5) 


\D 

00 
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To simplify notation, let 


i 


f = a 


Ar - 22. = 


A 




(C-6) 


Employing equation (C-6) and simplifying, quantity (C-5) reduces to 


fAR^'Jsin ji cos a - cos \i sin ajl^ 


(rtA) 2 R^r | + 2 [sin p cos n - cos p sin (ij 

f ' ' 

vhich can be further reduced by use of double angle quadrant relations 

* r ' P i 

and simplification, to 


b 

\ 2 t 


(rtA^R-jr 


^2 sin (p - jj.) + AR^ sin (pi - a)} ! (C-7) 


-Re substitution of equations (C-6) now provides the desirable form of 
term (C-7) after simplifying by angle quadrant relations, or 


(jtA) 2 %r 


{ 2 sin A(% - r) - AR^ cos A(R^ - r)^‘ 


(c-8) 


Thus, equation (C-2) may now be written as 


M™ = 2rt?5 / 

Fn Jr 


«o 


(jtA) nqr 


^2 ?in A(Ri - r) - ARj^ cos A(Ri - r)} 


J{_ 

\2r 


.W*± R c 


{ 2 . sin A(Rp - Rq) - ARp cos A^R^ - Rq)}- 


dr 
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Expanding the integrand of equation (C-9)* 'simplifying, and integrating 



After substitution of the limits and simplifying, the final relation 
for the generalized mass of the fuel is 



(C-12.) 


It is readily seen that expression (C-12) lessens considerably the 
labor involved in computing the effective inertia as opposed to the use 
of Bessel's functions. The eigenvalue problem for the approximate 
solution is developed as follows: 
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Substituting the equations (C-3.) and (C-4) into equation (2-24) 


for the appropriate order Bessel's functions yield 



(C-13) 


or, simplifying, (C-13) and rearranging yields 

AA jcos^ARo - sin (ar ± - |) - sin ^ cos (aRj. - 

+ 2XR ± |sin (aR q - 22.) cos ^7% - |) - cos ^AR 0 - sin " f) 

i 

+ 2AR 0 jsin (aRq - |) cos^ - 22.) - cos .(aRq - 0 sin (A% ^)j 

+ 4.fsin ^AR^ - ^) cos ^ARq - qp) - cos ^ARj_ - 0) sin .^ARq - ^— ) 

1 ' ‘ ' , / (c-14) 

Again, employing double angle relations, equation (-C-l4) simplifies to ^ 

A 2 R 0 R i sin A(Rq - Rq) + 2 ARq sin (aR q - ARj_ |) 

f ¥ 

+ 2 ARq sin (AR 0 - AR ± + |) + 4 sin A(Rq- R q ) = 0 


= 0 


or 


[A^Rq+^l sin A(R i -R 0 ) - 2AR£ cos A(Rq - Rq) + 2ARq cos ACRq - Rq) * 0 



A^R Q RjL + 4 


sin A(% - Rq) + 2A(R q - R ± ) cos A(.Ri - R Q ) = 0 (C-15) 


or, finally 


tan A(R^ - Rq) + 


2A(Rp - Rj_) 
k + R 0 R i ?v 2 


(C-l 6) 


which is the transcendental relation; developed from the approximations 
for the Bessel's function; for the roots A. 

23ie corresponding normalized mode shape is 
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Appendix D - Method of Obtaining Higher Modes 

In order to obtain modes -higher' than the fundamental in the 
solution of equation (3-48) an orthogonality sweeping procedure can be 
used and is developed as follows: 

Writing equation (3-48) for the gth and hth modes and rearranging. 


“g 2 t M ] ( y g ) - W ( y g) 

(D-l) 

W) - M <*.) 

(D-2) 


Premultiplying (D-l) by y^ and (D-2) by y g and subtracting the 
second equation from the first 

“ g 2 L y hJ M ( y g) - NM ( y g> 

■ ~ <% 2 L y gjl M ] ( y h) + L y sJt?] ( y h) = '° , 

* 

'"i 

By virtue of the symmetry of [m] and jjf], 

» i 

NM fe) = NW ( y &) . 

and a like relation, holds for [m] inserted for [if], thus permitting 
equation (D-5) to be simplified and rearranged to 

(cog 2 - <% 2 )[y g J[M] (y h ) = 0 (d-4) 

.and, for distinct . o>, then 

L y gjW ( y h} ® °> s ^ h Q* E *D* ( D -5) 



Thus, we have shown that any two eigenvectors of the system are orthogonal with respect to the 
inertia matrix M. To obtain the explicit relation required, relation (D~5) can written in 


expanded form as 






By performing, the indicated matrix multiplication } equation (D-6) takes the form 


H 

O 

-r 
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^0 


U)l F 0 + - + ^ h> + + ^Wo^Sj® 


+...+ |(<=i (g ' ) % i + J>2 (s) Bol» s 5o(s))]oi (h) 

+...+ ^Go^^o^s^oC 3 ) + ’ * * +< ^i^^o^ ::i Ps^i(s)j^ = 0 

Therefore, the solution form for any one of the components of the hth 
/ a o\^) 

eigenvector, say I 1 is 

a^(g) a^( +....+ a. n (g) an (^)x v . _ 


fo (h) 

R„ 


Rn Rf 


t f0 

F oR 0 


“~^y \ + Ag^^RohPg^CsJcQ^ +...+ 7v 2 (g)R 0 hp s C i (s)c i ^ h ^ 
c^S^RohPgCoCs)^^ +...+ Cx^Rohp^xCs^ 11 * 


or in more general terras 

y Q ( h ) = J=- U 1 yi (g)y i (h) + *** + 


(D-T) 


(D-8) 


Row utilizing the fact that- 


fjl 

R 0 


(*> a , (h) a Jh) 


Ro 


a n' 


Hr 


fn 

Ro 


(h> 


00 _ „ (n) 


= c n ' J , etc. 


(D-9> 


■which can be stated in general form as 


^(h) = y .(h) 
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A general expression can be written as’ iolldws 


or 


y 0 

(h) 

r 

0 

A l y l (g) 

yi 


0 

1 

• 

> - 

• 

• 

• 


• 

• 

■ 


• 

• 

1 


0 

- 

0 


r ^ 
■ y o 

y i 


(*)' 






{^} (h) - m w 


(D-10) 


where jjE>] is commonly termed the sweeping matrix and is applied as 
follows: From equation (3-5*0 the eigenvalue problem is 


therefore 


[H]{y} = A* (y) 

[H][s]{y> (h) =-V(y>^ 


(D-ll) 


and since the contribution of the first mode has been removed, iteration 
on equation (D-ll) will force convergence on the ' second mode and its 
corresponding frequency. For higher modes, the process is repeated, 
and the order of the sweeping matrix reduces for each successive 
operation. Note, however, that because of cumulative errors which are 
usually incurred in the successive sweeping matrix operation, it is 
wise to check the highest mode by inversion of the original matrix 
equation. - 



Appendix E - Computer Programs 


Those computer programs necessary to the solution of the systems 
problem developed in the investigation, contained herein and' developed 
specifically to support this investigation are presented on the pages 
following . 
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Propellant Cylinder Shear Modes 


PROGRAM A216S ( INPUT. OUTPUT* TAPE5= INPUT * TAPE6=OUTPUT ) 

DIMENSION ROOT ( 50 )»R(700)*PI700) * PPR ( 700) ♦XLAMBDA ( 1 00 > • R I J < 700 ) * 
1 R J ( TOO ) ♦ RBRJ ( 700 ) ♦ SUM ( 7 00.) * RY f 1 000 > * R I y { 1 000.) * T I TLE ( 8 ) . RBRY ( 700 ) 
COMMON RZ«Rl IC »AI « A2 ♦ A3* A4 * A5 * A6 * A7« AS ♦ 

real mfn«ifn« . r t 

EXTERNAL FOFX. 

NAMEL I ST/'NAM 1 XA « B « DX * RZ *R I « E l y E2 ♦ MAX-t * RBR , LMAX * RL *.DR « RMAX * GF ♦ 
IDELTA.EL*RC,FK*A1 *A2*A3*A4«A5*A6*A7*A8*GR* 

CALL CALCOMP, 

CALL LEROY* 

70 READ (5. 106 5TITLE* 

READ (5* NAM I ) , 

WR1TEC6*NAM1 ) , 

WRITE (6*1 05 ) « 

WR I TE C 6.* 1 06.) TI TLE ♦ 

SQ=SQRT(GF*GR/DELTA ) . 

P I =3« 14 159 * 

F A =2 • *P I *DELT A *EL /GR . 

IC = 1 

DO 30 I=1«LMAX* 

DDX=DX 

CALL ITR2«X.A,B,DDX*FOFX*Eri *E~2*MAXI * ICODE) « 

IF < I CODE ) 1 « 2 « I * 

2 XLAM0DA ( !• ) =X« 

A=X+OX* 

30 CONTINUE* 
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GO Tff 4. 

1 PRINT 1 02 * I CODE . I * 

102 FORMAT C1X.8HERROR N0,I4.28H IN ITR2 SUBROUT I NE * ROOT N0«I4>* 

IF 1 I *GT • I ).G0 TO 3 « 

STOP* 

3 LMAX=I-1 * 

4 DO 40 iFl.LMAX. 

NRR=0* 

K= 1 ♦ 

R<K>=RL* 

XLRO-XLAMBDA ( I ) *R I 
XLRBR=XLAMBDA ( I ) *RBR , 

. CALL BSSLStXLRO.RI J*2. IERR) « 

ME=1 * 

IF't I ERR ) 5 * 6 * 5 « 

5 PRINT lOO.ME.XLAMBOAC l ) ♦ 

LOO FORMAT (IX* 47HERROR IN’ B SSLS SUBROUTINE USING R ANO LAMBDA. AS . 13 « 
1E1S.7)* 

STOP. 

6 CALL BSSLS (XLRBR*RBRJ*2 * IERR) « • 

ME=3. 

•I F ( I ERR V5 « 8 * 5 * 

8 CALL BF4F(XLRO.RlY*2, IERR.-.1 1. 

M3=l. 

IF( IERR.) 9,:l 0.9. 

9 PRINT 10L* I ERR * M3 * XLAMBDAI I ) * 

10 1 FORMAT Cl X*8HERR0R NO.T4.42H IN BF4F SUBROUTINE USING- R AND LAMBDA 
1 AS. I3.E15.7) , 

STOP. 

10 CALL BF4F(XLRBR-.RBRY,2.-IERR.'-1 > « 

M3=3. 

IF( IERRJ9..12.9. 

- 12 Ff=A3*XLRO*R7J( 1. ) -A4*R,I J (2 >... 

F2=F1*RBRY<2 ) « 

F3=A7*XLRO*RIYtl )-A8*RIY<2>, 

F4=F2-F3*RBRJ ( 2 ) . 

18 XLRs=XLAMBDA (I )*R(K) . 

.WR I TE'( 6 ♦ 1 1 0 ) XLAmBDA ( I ) * R « ) *XLR . 

110 FORMATC* lambda =*EL5.7.* R =*E15.7.* XLR =*E15.7>« 

CALL BSSLS 1XLR .RU-. 2 * IERR ) * 

ME=2 ♦ 

I F ( I ERR ) 5 » 1'3 . 5 . • 

13 CALL BF4F CXLR.RY *2* IERR * — 1 ) » 

M3=2. 
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IF 1 1 ERR )9»14*9« 

14 P!K)=(F1*RY!2)~F3*RJ!2) >/F4, 

RDR=R (K )+DR* 

IF-CNRR.EQ.l >60 TO 20, 

IF'!RDR,LT.RMAX> 15,50, 

15 K=K+1. 

R!K>=RDR, 

GO TO 18, 

50 NRR=1, 

K=K+1, 

R(I<C,)=RMAX, 

GO TO 18, 

20 DO 1'6 M=1,K, 

16 PPR!M)=R!M)*P!M)**2, 

CALL TRAP (PPR,R,K,SUM»0 ) , 

MFN=FA*SUM CK ) , 

OMEG=XL.AMBDA 1 1 ) *SQ 
EF=OMEG/!2.-*PI > , 

•IFN=RZ**2*MFN . 

BUF= I FN*OMEG**2 , 

PROD=RC*P!K> 

WR I TE (6,1 07 ) , 

WRITEI6, 108JXLAMBDA ( I > , OMEG , EF , MFN , I FN , BUF , PROD , 

CALL PLOT ! K * R » P , EF , FK ) , 

40 CONTINUE* 

GO TO 70, 

105 FORMAT! 1H1 ,/15X«*ELASTIC FUEL SHEAR TWIST MODES*) 

106 FORMAT ( 8A 10) , 

107 FORMAT !/V *5X ,6Hl_AMBDA * 9X»SHFREQ RPS *7X, 8HFREQ CPS»7X*8HGEN MASS, 
1 1 OX , * I NF* ,1- 1 X , *BUF* ,-i OX , *PROD* > « 

108 FORMAT! 7E1 5,7 ) , 

END, 

FUNCTION FOFXCX) 

DIMENSION RANS! 700), RIANS! 700 ) ,RST !7O0 > » RI ST! 700 ) , 

COMMON RZ,RI * IC,A1 ,A2,A3,A4 ,A5,A6, A7,A8, 

XR=X*RZ» 

XR I =X*R ir, 

CALL BSSLS < XR, PANS *2,^1 ERR )’» 

I.F ( I ERR ) 1 ,2,1 , 

1 PR I NT 10,1.x, 

101 FORMAT! IX, 36HERROR IN BSSLS SUBROUTINE, FOFX ROOT ,E1 5 ,7 ) 

STOP, 

2 CALL BSSLS! XR I ♦ R I ANS ,2,1 ERR ) , 

IF! IERR) 1 ,3,1 , 
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3 CALL BF4F(XR*RST*2* IERR'*-1 ) ♦ 

IF ( I ERR )5i4t5« 

5 PRINT 1 02 * I ERR * X ♦ 

1*02' FORMAT CIX*8HERR0R NO, I 4 *311-1 IN BF4F SUBROUTINE* FOFX ROOT,E15#7>* 
STOP* 

4 CALL BF4FCXRI *RIST*E, IERR*-1 ) * 

IF< IERR)5.*6*5* 

6 FOFX= CA1*XR*RANSU )-A2*RANS(2 ) ) * ( A3*XR I *R I ST C 1 )-A4*RISTC2> )-CA5*XR 
1/*RST ( 1 > -A6*RST < 2 > ) * < A7*XR I *R TANS ( 1 > - A8*R I ANS ( 2 ) ), * 

I F-( I C*EQ • 1 ) 20 , 30t* 

20 WR I TE (6*103)X*F OFX ♦ R I ANS < 1 >*RI ST-C 1 ) * 

103 FORMAT (yiX.*2HX= *Ei5.7<*2X*3HFX=*El’5*7*2X*7HJRI Cl )=*E15»7* 7HYR I < 1 ) = 
1E15.7/) 

, IC=IC+1 • 

30 CONTINUE* 

- RETURN* 

END ♦ 

SUBROUTINE TRAP (Y *X *NL * SUM « SUM1 )* 

DIMENSION Y ( NL ) * X < NL ) * SUM ( NL ) * 

SUM ( 1 ) =SUM 1 * 

NLM=NL— 1 * 

DO .1 N= 1*NLM* 

NP=N+1 * 

1 SUM(NP)=SUM(N)+CYfNP)+Y(N) >* < XTNP )-X CN > )/2. ♦ 

RETURN* 

END 

SUBROUTINE PLOT (K* R *P* EF* FK ) * . 

DIMENSION R(700.) *PC700) **(702) *Y (702 V* PH I C700) * 

DO 10 I -1 *K* 

10 PHI C I )=PU )/RC I )* 

PMAX=PH1 (1 ) 

DO 20 1=1 *K ♦ 

IF (PHI ( 1 >,• GT *PMAXM 5* 20 ♦ 

15 PMAX=PH 1 (I) * 

20 CONTINUE* 

F=1 •/<FK*PMAX) * 

DO 40 1=1 *K* 

PHI (IJ=PHirt >;*F* 

X( I)=R( 1 )*COSCPHl (1 ) ) * 

40 YC I )=RC I )*SINCPHI < I M * 

WRITE ( 6*1 00 ) (RT-t hPCI ) *PH1 ( I ) * X H 1 * Y ( I ) ♦ I “ 1 * K.) * 

1 00 FORMAT < //7X* *R* * 1 4X * *P* * 1 2 X * *PH I * ♦ 1 4X * *X* * 1 4X * *Y*/ C 5E 1 5 • 7 > ) ♦ 

CALL GRIDCO* *0* *1, •* 1* *5*5).* 

CALL ASCALE CX*5* *K* 1 ,20# ) * 
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CALL ASCALE (Y.5..K, 1 » 20 » ) « 

’ K1 =K+1 ♦ 

' K2=K+2, 

YCK2 )=XCK2) ♦ 

CALL AXES10. ,0. ,C.»5. .X(K1 ) ,X(K2 > .1 .-,2. . 1HX, • 125,-1 > « 

CALL AXES<0. ,0. ,90. » 5. , Y (Kl ) ,-YCK2) » 1 . ,2. » 1HY. • 125.1 1 . 

CALL N0TATEC2. ,5«5«. 14^11 HFREQUENCY =,0,,11>» 

CALL NUMBER (3. 8, 5. 5.. 14, EF,0..3), 

CALL -LI NPLT < X, Y,K »-l , 0, 0 ,0, 0 ) , 

CALL CALPLT (10*«0.« —3 ) i 
RETURN , 

END, 

PROGRAM A2168 CASE 7CK-1A FREE-FREE 4-29-70 

SNAMl A=.0I ,8=20. ,DX=«01 , RZ=7 • ,R I = ,'5 ,E1'= .1 E— 6, E2= . 1 E— 6»MAx I ®100,RBR=7 • , 
GR— 386. , . 

LMAX=15,RL=.5,DR=.25,RMAX=7, ,GF=90.. ♦ DELTA=53.460E--3,EL=1 • ,RC=.5,A1=1 * ,, 
A2=2. , A3= 1 . , A4=2 . , A5= 1 * , A6 = 2 . , A7.= l , , A8=2. ,FK~5 .73S 


— * 
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PROGRAM AO 107 ( INPUT « OUTPUT ♦ TAPE5= I NPUT, TAPE6=0UTPUT ) 

DIMENSION XM (40,40 > «W( 40 ) ,XN (40,40 ) ,X IF (40 ) ,PP (40 ) ,V(40 ) »ZI (40 ) , 
1XG(40),WI(40).P(40>« 

2 , FK ( 40 ) *WV ( 40 ) ,GFV(40> ,R(10)» 

DIMENSION HEADO). 

- NAMELIST /NAMI /W»XL ,GF * RZ* PP*P *V»WI«ZI ,XH,RH0,NN, I I , X IF« FI ,14, 
III, 16, TEST, G, Cl ,ER I ,XG , RN, FK , TA3 , DTA3 ,TF , NROOTS , EPS » ER2 , 

1 READ <5, 103 )HEAD , 

IF (EOF ,512,3, 

2 STOP, 

3 READ (5^ NAMI ) , 

WR 1 TE ( 6 * NAMI ) , 

WRITE(6« 104 1HEAD, 

DO 10 1=1,40, 

DO 10 J=1 ,40, 

XM(I«J>=0., 

10 XN t 1 « J)=0, » 

Pt=3. 14159, 

NAB=NN+I 1+2, 

NC=NN+II+1 , ■ 

00 50 J=l ,NN, 

IK(W(J> .EQ*0. >22,24, 

22 R( 1 >=0. » 

NR=1, 

GO TO IS, 

24 B=-FK(U)*W( J)**2/GF« 

C=-W(J>**2, 

'K=3. 

T3=TA3« 

DO IS 1=1, NROOTS* 

8 CALL SEARCH2(T3,DTA3,TF,ER2.FN2,K,DFN2) , 

GO T0(4»5,6),K, I 

4 FN2=T3**2+B*T3**EPS+C, 

DFN2=2*T3+EPS*B*T3**(EPS-1 > , 

GO TO 8, 

5 R( I 1 =T3 , 

IS CONTINUE, 

NR=NROOTS, 

,GO TO 18, 
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6 NR= I— 1 « 

IF(NR.NE.O>GO TO 18, 

WR I TE ( 6 » 1 06 ) W ( J ) ♦, f 

106 FORMAT t/* NO ROOTS FOUND .FOR. W ' =*E1 5*7 > ♦ ' 

stop*. ■ ; ! 

18 WR1TE(6,105)W( J ) * (R ( 1) * 1 = 1 
105 FORMATS /* VI =*E 15*7* /* ROOTS ‘ =*.( 6E15.7M * . 

WV ( J ) =R 111* 

I F ( EPS • EQ * 0 » >26.28* 

26 FA = 1 .* 

GO TO 32, 

28 FA=R(1)**EPS. 

32 GFV(J-)=FK(J)*FA, 

WRITE (6.1 07>GFV{ J> . 

107 FORMAT!.* GFV =*E15.7)» 

50 CONTINUE. 

DO 20 N=1 ,NN . 

- XM (N, N)=X IF ( N) . 

XN(N,N)=WV(N>**2*XIF(N)+F1*2*PI*GFV(N)*(RZ**2)*PP<N>**2. 
XN ( N ♦ NC )'~RZ*P < N.) , 

XN(NC.N-) =RZ*PCN> , 

20 XN t NAB , N ) =XN (N , NAB ) =RZ*GFV ( N > * < PP ( N > -P ( N > /RN > , 

NCI =NN+ 1 1 , 

NN1=NN+1'«. 

1 = 1 , 

DO 25 N=NN1,NC1, 

XM(N.N)=V(I). 

XNIN,N)=V(I)*WI (t )**2. 

XN ( N . NC > =— RN*Z 1(1). 

XN(NC, N ) =—RN*ZI (T 1 , 

XM ( N « NC+ 1 ) =XM ( NC+ 1 , N > =RZ*XH*RHO*Z I ( I )•, 

25 I =1+1, 

DO 30 JJ=1 ,NA8 , 

WRITE (6. 100 >JJ, (XM( JJ. J ) , J=1 .NAB ) . 

30 CONTINUE. 

WRITE(6,101 ) , 

DO 40' JJ=1 ,NAB . 

WR I TE (6,102) J J . ( XN ( J J ♦ J ) , U= 1 , NAB > » • 

40 CONTINUE, 

CALL. INVERT (NAB, 14, II , 16.40 .TEST ,G»C:1 ,ER1 .XM.XN.XG.) , 

GO TO 1 ♦ 

100 FORMAT ( * M(*,I2,.*,1 >.=*/ (8E1 S.7 > V* 

101 FORMAT (//) , 

102 FORMAT (* >NC*»I2,*, 1 ) =*/ CBE1 5.7 > > ,. 
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1-03 FORMAT ( 8A 1 0, > * 

104 FORMAT (/* 1 0X*3A 1 0 ) * 

END. 

SUBROUT I NESEARCH2 ( TA2 « DTA2 *TF « ER2 ♦ FN2 « K , DFN2 >.* 
REALLAMA,LAMB', 

GOTO (t«2»3)*K* 

1. 1F(JF.GE.2)G0T04«- 

5 TA1 =TA2* 

FN1 -FN2* 

DFN1-DFN2, 

IE=0« 

FS=0., • 

GOT06* 

2 FNI = (DFK}2/'A0S‘(OFN2 ) )•*• 00001 * 

7 DFN1=DFN2. 

TA1^TA2«. 

JF=1 t 
G0T08* 

3 DFN1=0*« 

FNt=0. t 
TAI=0. « 

JEND=0. 

JF=0 ♦ 

8 LAMA=0.« 

LAMB=Oo ♦ 

JW=0« 

JK=0 * 

6 T A2=T A2+DT-A2 * 

JF= JF+ 1 « 

9 I F t TA2*GT «TF JGOTOl 0* 

11 K=1 * 

RETURN* 

10 IF< JEND*EG*3.)G0T012« 

TA2=TF« 

JEND= JEND+ 1 
GOTO 11 * 

12 K=3 ♦ 

RETURN * 

4 IFCJK.EQ.l JG0T013* 

IF (FN2/FN1 ) 1 4* 15*16* 

16 IFIIE*EQ*0 ) GOTO 17* 

18 I F'( (FN2—FS >/FN2*LE*ER2 IGOTOl 9 • 

A1=DFN1 i 

A2= (FN2-FN1 — DFN,1* (TA2— TA1 ) )/<TA2-TAl >**2« 



116 


TA3=TAl-AlV(2.*A2) • 

FS-FN2» 

IFCIE.EQ.l 1GOT020. 

IE=1 « 

G0TO25. 

20 TA1=TA2« 

FN1=FN2« 

DFN1.=DFN2 * 

25 TA2=TA3. 

G0T09. 

19 FN2=SFN2* 

•DFN2=SDFM2. 

TA2=STA2» 

FN1=FN2. 

IE=0. 

FS=0 • « 

G0T07. 

I'7 I F ( OFN 1 /DFN2 • GE • 0 o I.GOT05 • 

! F t DFN 1 *FN 1 . GE. 0 • > G0T05 ♦ 

SFN2=FN2» 

SDFN2=DFN2 1. 

STA2=TA2« 

. G0T018» 

14 IE=0, 

‘FS=0 • ♦ 

13 A0=FN1 « 

•A 1 =DFN1 . 

A2= (FN2— FN1 *rDFNt * ( T A2— TA1 ) >/(TA2-TAl >**2» 

T A3 1 = ( -A 1 +SQRT ( A 1 **2-4 . *A0* A2 > ) / C 2 . *A2 M-TA 1 
TA32=(-AI-SQRT(-A1**2-4.*A0*A2) >/(2.*A2>+TAl « 
IF ( T A t « LE • T A3 1 .AND.TA31 .LE.TA2 JGOT022. 

IF (TA1 ,LE.TA32.AND.TA32.LE.TA2 >G0T023« 

Q=ABS 1TA31 -TA2 » , 

R=ABS CTA32-TA2 ) . 

T.A3=AMIN1 CQ.R) . 

TF (TA3.EQ.Q 1GOT022* 

G0T023. 

22 TA3=TA31'. 

G0T024. 

23 TA3-TA32* 

24 LAMA=LAMB « 

LAMB=TA3. 

JK=1.» 

JW=JW+1 » 
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IF < JW.LT.2>GOTO20-. 

TEST='< LAMB-LAMA )/LAMB, 

IF TABS (TEST ) .GE.ER2 IG0TO20 » 

- S3 K=S., 

TA2=LAM0, 

RETURN, 

15 LAMB=TA2. 

GOT 021 ♦ . , , ’ . 

• 102 FORMAT <1X* INTERPOLATED TIME IS NOT BETWEEN ,T1 AND ,T2* > i* 

END. ' . - 

SUBROUTINE INVERTtN. 14, Il», 1 6 .NMAX.TEST ».G, Cl.'.ERl . A , B. XG )-, 

C PROG. NO. 0107- RCM ROSS GOBLE. , ' > V. 

EX-TERNALMAT I NV « ’if. ’ * * „ 1 

C IL=LENGTH OF GUFSS VECTOR. t * , 

C 11=0 NORMALIZE ON MAXIMUM ELEMENT ♦ =N0 »> NORMAL I ZE- ON NO., 

C I 2= ELEMENT WHICH IS NORMALIZED ON. 

C I4=N0. OF MODES WANTED-. ■ * 

DIMENSIONA <40,401 ,XB <40 > ,X< 40 > ♦ AH <40 ,40 ) ♦ BC40, 40 ) »C <40.40 ) , 

1 INDEX.<40 ,2 > .AM < 40, 40 > , XV<40 , 40 ) ,WB<40.4 > « II <31 .ERR (4 ) , 

1 XG <40 ) » 

II <1>=I1 
II <2>=I4 
I I < 3.) = 1 6 
ERR < 1 V=TEST 
ERR' <2 )=ERl 
ERR <3 )=C1 
ERR <4 >=G- 

14 D01SI=1»N» 

DOISJ=l,N-, 

15 AM < I ♦ J > =A <-I ♦ J I , 

' CALLEIGEN< I I .N.NMAX. 13. ERR ♦ SUM 1 , A , B.AM .AH.C, X-V, INDEX ♦ XG » X ♦ XB .'WB , 

1MATINV) 

WR I TE (6.110) SUM 1*. 

WRITE <6,1X2 ) • 

‘ WRITE<6.1 13) . 

D01 1=1 .13, 

WRITE<6.,109) . 

WR ITE <6 ♦ 111 ) <WB<I,J>,J=1«4> 

. WRITE<6.109) ♦ 

1 WR I TEX 6 « 1 06 ) < XV < J , I ) » J= 1 » N ) ♦ 

RETURN. 

100 FORMAT < 1X1 1A6) « 

106 FORMAT < 7E1 6. 8) « 

109 FORMATC IX/ ) « 
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no FORMAT < 1X4HSUM=.E 16*8) « 

111 FORMAT (4E16. 8). 

1 T2 FORMAT ( 6X3HLAM 14X1 HWi 5X3HSUM1 0X2HI 5 ) ♦ 

113 FORMAT { 7X4HX.( !'»■)•« 

END. 

SUBROUT I NEE1 GEN < I I . ! L « NX . L3. ERR .SUM 1 . A.B .MASS » AH .C.XV. INDEX'.XG.X. 
1XB.WB.MATINV) 

REALMASS . LAM . LAmX 

D I MENS 1 ON A (NX. It.) «B(NX. IL) .MASS (NX. IL).AH(NX. IL ) .XV(NX . 1L) .XG( ID 
1 .X( I L ) « XB ( I L ) » WB ( NX * 4 ) « C ( NX ♦ I L ) . I NDEX’ ( NX . 2 ) ♦ I I (3) .ERR (4) 

11=11(1) 

14=11(2) 

I 6=1 1(3) 

ER=ERR ( 1 ) 

ER1 =ERR (2 ) 

C1=ERR(3) 

G=ERR(4> 

D04 1=1 . IL. 

D04J=1 . IL 

4 B ( I * J >=8 ( I « J )+Cl *MASS ( I ♦ J > ♦ 

DO 60 1=1 .IL. 

WRITE(6.100)I. (B(I.J).J-l. ID. 

60 CONTINUE. 

100 FORMAT (* ADD (*• * 12 .* «T ) =-*/ (8ET5 » 7 ) ) . 

DOl 1 (=1 . IL. 

DO IT J=1 , IL. 

11 B(T,J) = B(I.J-)/G, 

CALLMAT I NV ( B . I L ♦ C . O’. DETERM . X . I NDEX . NX «T SCAL E ) 

D022I=1 . IL* 

D022K=.l . IL. 

C ( I «K) = 0 ... 

D022 J=1 ♦ IL. 

22 C ( I .K)=B(!.J)*MASS'TJ.K)+C(I.K)» 

SUM 1=0. * 

00231=1 . IL. 

D023K=1 . IL. 

23 SUM1~SUM1+C(I.K)*C'(K.I ) . 

IF (Cl *EQ«0* )G0T071 . 

TEST1 =ABS (G/Cl ) . 

71 13=0 . 

SUM=0. . 

LAMX=0.« 

DO20 1= 1 « IL . 

•D020J=1. IL-. 
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20 B ri«J)*0" 

0021 I -1 * I L* 

J= It 

21 B 

00191=1 * I L * 

D01>9J=1 4 IL* 

19 AH C I * J ) =C ( 1 * J) « 

26 rs=o* 

D031 1=1 *IL* 

31i XB(I)=XGU>* 

8 002I=l4lL4 
, X< 1 >=0. 4 
D02J=1 * IL* 

2 X(I )=AH.<I*J)*XB<J>+XCI > ♦ 

15=15+1 4 

I F ( 1,5* LT • ( 16-10 )->G0T050t- 
WR I TE ( 6 4 1 03 )LAM f 
50 TF( I5*UE. I6X50T0254 
WR I TE ( 64 1'05 ) • 

RETURN* 

25 I-F.(T 1 •NE#0)G0T01 * 

12=1 4 

0031=1 « IL * 

IF ( ABS (X ( 12$ ).*>GF • ABS ( X ( I ) > 1G0T03* 
1 2= I * r 

3 CONTINUE * 

G0T05 • 

1 12=11 * 

5 LAM=ABSOOC I2 >> * 

D06I=1 ♦ IL* 

6 XtJ >=XC I J/LAM4 
DO 181=1 4 IL* 

18 XB ( t ) =X C I )-* 

TLAM=ABS t ( LAM-LAMX > /LAM ) * 
IF<TLAM*LE#ER>G0T07* 

,LAMX=LAM* 

G0T08 * 

7 Q=ABS( (G/LAM)-Cl ) * 

SUM=SUM+.LAM**2* 

W=SGRT ( Q ) « 

rFICl • EQ • 0 * ) GOT070 * 

■IF CABS ( LAM-TEST 1 )*GT*ER1 >G0TO70* 

w=o* « 

70 13=13+1 • 
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WB C 13* 1 )=LAM 
W8C I3*2)=W 
W8C I3.3)=SUM 
WB( 13*41=15 
D09I = 1.« IL* 

9 XV ( I « 1 3 ) =X t I )♦ 

IFC I3.LT. I4.)GOTO30. 

RETURN * 

30 00101=1 • IL • 

D01 0K=1,.« IL* 

A C I «K ) =0 • « 

DOl 0J=1 « IL.* 

tO Ad *K)=XCI >*XCJ)*MASSCJ.K)+AC.I*K>« 
00121 = 1 «,IL* 

XB CI>=0.« 

DOI2J=l « IL* 

1'2 XB C I)=XCJ>*MASS'{J*n+XB CI>* 
DEN0=XB <1>*X<1')« 

D013I =2 ♦ I L . 

1*3 DEN0=DEN0+X8 C1)*XCI). 

00141=1 * IL* 

001'4J=.| ♦ IL « 

14 ACI*J)=ACI * JJ/DENO* 

00151 = 1 * IL* 

DO!5J=l ♦ IL* 

15 A < I * J) =B < I « J )LA 1 1 » J ) « 

00161=1 *IL* 

0Ol'6J=*l ♦ IL* 

1*6 B <I.J)=A.(I*J>« 

DOl 71= 1 * IL*. 

0017K=1 * IL « 

AH (-1 * K )=0»’* 

OOl 7J=1 *>IL 

17 AH Cl *K)=Ccr* J>*ACJ»K)+AHC I *K) * 
GOT026 « 

1 03 FORMAT C 7E 1 6 « 8 > * 

1*05. FORMAT Cl XI 7HD0ES NOT CONVERGE)* 
END* 



